#Replication material for the Generalizability of IR experiments beyond the U.S.

##Load Relevant Packages----
#make sure latest version of R is installed
rm(list=ls())
library("tidyverse")
library("broom")
library("estimatr")
library("ggplot2")
library("metafor")
library("ggthemes")
library("ggh4x")
library("haven")
library("ggimage")
library("devtools")
install_github("jimjam-slam/ggflags", dependencies = TRUE)
library("ggflags")
library("texreg")
library("ltm")
library("interflex")
library("stargazer")
library("tidytext")
library("wordcloud")
library("gridExtra")
library("quanteda")
install_github("naoki-egami/exr", dependencies = TRUE)
install_github("naoki-egami/evalid", dependencies = TRUE)
library("evalid")
library("exr")
library("hettx")
source("simulations/functions.R")

set.seed(1)

## Load Simulated data-----
###Read in data----
#USA
usa_data <- read.csv("data/survey_data/usa.csv")
usa_data <- usa_data[3:nrow(usa_data),]
usa_data <- as.data.frame(usa_data) %>% mutate(., country = "USA")

#India --- 
india_data <- read.csv("data/survey_data/india.csv")
india_data <- india_data[3:nrow(india_data),]
india_data <- as.data.frame(india_data) %>% mutate(., country = "India")

#Germany
germany_data <- read.csv("data/survey_data/germany.csv")
germany_data <- germany_data[3:nrow(germany_data),]
germany_data <- as.data.frame(germany_data) %>% mutate(., country = "Germany")

#Brazil
brazil_data <- read.csv("data/survey_data/brazil.csv")
brazil_data <- brazil_data[3:nrow(brazil_data),]
brazil_data <- as.data.frame(brazil_data)%>% mutate(., country = "Brazil")

#Israel
israel_data <- read.csv("data/survey_data/israel.csv")
israel_data <- israel_data[3:nrow(israel_data),]
israel_data <- as.data.frame(israel_data)%>% mutate(., country = "Israel")

#Japan
japan_data <- read.csv("data/survey_data/japan.csv")
japan_data <- japan_data[3:nrow(japan_data),]
japan_data <- as.data.frame(japan_data)%>% mutate(., country = "Japan")

#Nigeria
nigeria_data <- read.csv("data/survey_data/nigeria.csv")
nigeria_data <- nigeria_data[3:nrow(nigeria_data),]
nigeria_data <- as.data.frame(nigeria_data)%>% mutate(., country = "Nigeria")


all_countries <- bind_rows(usa_data, india_data, germany_data, 
                           brazil_data, israel_data, japan_data,
                           nigeria_data)


###Recode variables-----

all_countries$AGE<-ifelse(all_countries$AGE=="",NA,all_countries$AGE)
all_countries$S<-ifelse(all_countries$S=="",NA,all_countries$S)
all_countries$ideology<-ifelse(all_countries$ideology=="",NA,all_countries$ideology)
all_countries$education<-ifelse(all_countries$education=="",NA,all_countries$education)
all_countries$voting_notWVS<-ifelse(all_countries$voting_notWVS=="",NA,all_countries$voting_notWVS)

all_countries_clean <- 
  all_countries %>% 
  mutate(.,
         #Recode main treatment variables
         DP_treat = case_when(
           dem_condition == "democracy"~1,
           dem_condition == "non_democracy"~0),
         
         AC_treat = case_when(
           aud_condition == "stay_out"~0,
           aud_condition == "not_engage"~1),
         
         IL_treat = case_when(
           il_condition == "law"~1,
           il_condition == "no_law"~0),
         
         IPE_treat = case_when(
           ipe_condition == "easier"~0,
           ipe_condition == "harder"~1),
         
         #for additional AC analysis, we also break treatment by:
         belligerence = case_when(
           aud_condition == "stay_out"~0,
           aud_condition == "engage"~1),
         
         inconsistency = case_when(
           aud_condition == "engage"~0,
           aud_condition == "not_engage"~1),
         
         #Recode main outcome variables:
         
         #outcome1_dem: original - attack country's site, high # indicate less support
         #outcome1_dem_new: support attack, high # indicate more support
         # "By attacking the country's nuclear development sites now, Nigeria could prevent the country from producing any nuclear weapons. Would you favor or oppose using Nigeria's military to attack the country's nuclear development sites?"
         
         outcome1_dem_new = 6-as.numeric(outcome1_dem),
         
         #outcome2_dem: original - join operation, high # indicate less support
         #outcome2_dem_new: support joining military operation
         # "By joining a joint international military mission now, Nigeria could prevent the country from producing any nuclear weapons. Would you favor or oppose using Nigeria joining an international military mission?"
         
         outcome2_dem_new = 6-as.numeric(outcome2_dem),
         
         #outcome_aud: original - high # indicate less approval
         #outcome_aud_new: support leader (high # more support)
         # "How much would you approve or disapprove of the way [leader] handled the situation?"  
         outcome_aud_new = 8-as.numeric(outcome_aud),
         
         #outcome_il: original - oppose torture, high # less support for torture
         #outcome_il_new: original - support torture, high # more support for torture
         # "If things happened just as we described, to what extent would you agree or disagree with the following statement: [country] should use interrogation methods involving torture on captured combatants."
         outcome_il_new = 6-as.numeric(outcome_il),
         
         #outcome_ipe: original - support harder barriers, high # indicate support to make it harder to purchase [country's] companies
         #outcome_ipe_new: same as original
         # "Should [country] make it easier or harder for companies from this country to buy [country's] companies?"
         outcome_ipe_new = as.numeric(outcome_ipe),
         
         #Recode manipulation checks:
         
         #manipulation_dem: "1" is democracy, "2" non-dem
         DP_man = case_when(manipulation_dem=="1"~1,
                            manipulation_dem=="2"~0),
         
         #manipulation_aud: "1" is stay_out, "2" not_engage, "3" engage.
         AC_man = case_when(manipulation_aud=="1"~0,
                            manipulation_aud=="2"~1,
                            manipulation_aud=="3"~0),
         
         #manipulation_il: "1" is law, "2" is no law
         IL_man = case_when(manipulation_il=="1"~1,
                            manipulation_il=="2"~0),
         
         #manipulation_ipe: "1" is easier, "2" is harder
         IPE_man = case_when(manipulation_ipe=="1"~0,
                            manipulation_ipe=="2"~1),
         
         #Recode plausibility qs to make higher #s more plausible
         plausible_DP = 5-as.numeric(plausibility_dem),
         plausible_AC = 5-as.numeric(plausibility_aud),
         plausible_IL = 5-as.numeric(plausibility_il),
         plausible_IPE=5-as.numeric(plausibility_ipe),
         
         plausible_DP1 = case_when(
           plausibility_dem == "1"~ "Plausible",
           plausibility_dem == "2"~ "Plausible",
           plausibility_dem == "3"~ "Not Plausible",
           plausibility_dem == "4"~ "Not Plausible"),
         
         plausible_AC1 = case_when(
           plausibility_aud == "1"~ "Plausible",
           plausibility_aud == "2"~ "Plausible",
           plausibility_aud == "3"~ "Not Plausible",
           plausibility_aud == "4"~ "Not Plausible"),
         
         plausible_IL1 = case_when(
           plausibility_il == "1"~ "Plausible",
           plausibility_il == "2"~ "Plausible",
           plausibility_il == "3"~ "Not Plausible",
           plausibility_il == "4"~ "Not Plausible"),
         
         plausible_IPE1 = case_when(
           plausibility_ipe == "1"~ "Plausible",
           plausibility_ipe == "2"~ "Plausible",
           plausibility_ipe == "3"~ "Not Plausible",
           plausibility_ipe == "4"~ "Not Plausible"),
         
         #Recode info leak variables "1" ~ thought of country, "2 didnt
         DP_info = case_when(info_leak_dem == "1" ~ 1,
                             info_leak_dem == "2" ~ 0),
         AC_info = case_when(info_leak_aud == "1" ~ 1,
                             info_leak_aud == "2" ~ 0),
         IL_info = case_when(info_leak_il == "1" ~ 1,
                             info_leak_il == "2" ~ 0),
         IPE_info = case_when(info_leak_ipe == "1" ~ 1,
                             info_leak_ipe == "2" ~ 0),
         
         #Recode info leak variables "1" ~ country, "2" - no country
         DP_info1 = case_when(info_leak_dem == "1" ~ "Country",
                             info_leak_dem == "2" ~ "No Country"),
         AC_info1 = case_when(info_leak_aud == "1" ~ "Country",
                             info_leak_aud == "2" ~ "No Country"),
         IL_info1 = case_when(info_leak_il == "1" ~ "Country",
                             info_leak_il == "2" ~ "No Country"),
         IPE_info1 = case_when(info_leak_ipe == "1" ~ "Country",
                              info_leak_ipe == "2" ~ "No Country"),
         
         #Recode moderators
         #democratic norms - recode more support
         dem_norm1 = as.numeric(democratic_norms_1),
         dem_norm2 = as.numeric(democratic_norms_2),
         dem_norm3 = 6-as.numeric(democratic_norms_3),
         dem_norm4 = 6-as.numeric(democratic_norms_4),
         
         #hawkishness - recode more hawk
         hawk1 = as.numeric(hawkish_attention1_5),
         hawk2 = 6-as.numeric(hawkish_attention2_2),
         hawk3 = 6-as.numeric(hawkish_attention2_3),
         
         #legal obligation - recode more obligated
         legal1 = 6-as.numeric(legal_obligation_1),
         legal2 = 6-as.numeric(legal_obligation_2),
         legal3 = 6-as.numeric(legal_obligation_3),
         
         #create FEs
         country_fe = as.factor(country),
         gender_fe = as.factor(S),
         education_fe = as.factor(education),
         voting_fe = as.factor(voting_notWVS),
         age = as.numeric(AGE),
         #more covariates
         fear = (as.numeric(fear_1) + as.numeric(fear_2))/2,
         ideology = as.factor(ideology),
         #regroup education into university_level (nigeria had additional categories)
         university = case_when(
           country=="Nigeria"&education=="6"~1,#university with degree
           country=="Nigeria"&education=="10"~1,#MA
           country=="Nigeria"&education=="11"~1,#PhD
           country=="Nigeria"&education=="5"~0,
           country=="Nigeria"&education=="7"~0,
           country!="Nigeria"&education=="5"~1,#BA
           country!="Nigeria"&education=="6"~1,#MA
           country!="Nigeria"&education=="7"~1,#PhD
           education=="1"~0,#no highschool
           education=="2"~0,#some highschool
           education=="3"~0, #highschool grad
           education=="4"~0) #some
         )

##Scale all main outcomes for comparison
all_countries_clean<-
  all_countries_clean%>%
  group_by(country)%>%
  mutate(outcome1_dem_scaled = scale(outcome1_dem_new),
         outcome2_dem_scaled = scale(outcome2_dem_new),
         outcome_aud_scaled = scale(outcome_aud_new),
         outcome_il_scaled = scale(outcome_il_new),
         outcome_ipe_scaled = scale(outcome_ipe_new))%>%
  ungroup()


#index moderators
#dem_norms
all_countries_clean$dem_norm <- (all_countries_clean$dem_norm1 + 
                                   all_countries_clean$dem_norm2 +
                                   all_countries_clean$dem_norm3 +
                                   all_countries_clean$dem_norm4)/4

#hawk
all_countries_clean$hawk <- (all_countries_clean$hawk1 + 
                                   all_countries_clean$hawk2 +
                                   all_countries_clean$hawk3)/3

#legal
all_countries_clean$legal <- (all_countries_clean$legal1 + 
                               all_countries_clean$legal2 +
                               all_countries_clean$legal3)/3

#create variable for resp who took the survey in english
all_countries_clean$english<-ifelse(all_countries_clean$UserLanguage=="EN",1,0)

#create minutes var
all_countries_clean$min<- as.numeric(all_countries_clean$Duration..in.seconds.)/60

#Plot response latency----

summary_minutes <- all_countries_clean%>%filter(ideology!="")%>%
  group_by(country) %>% 
  summarise(mean.value = mean(min))

all_countries_clean%>%filter(ideology!="")%>%
  ggplot(., aes(x=min))+
  geom_density(alpha=0.4,fill="#FF6666")+
  facet_grid(~country)+
  geom_vline(data = summary_minutes, aes(xintercept = mean.value), colour = "red") +
  geom_text(data = summary_minutes,aes(label= paste("mean =", format(round(mean.value, digits = 1))), colour= "red", y=0.1,x=70)) +
  xlim(0,110)+
  xlab("Minutes to complete survey")+
  ylab("Density")+
  theme_bw()+
  theme(strip.background =element_rect(fill="dodgerblue4"),
        strip.text = element_text(colour = 'white'),
        panel.grid.major = element_blank(),
        legend.position = "none",
        strip.text.x = element_text(size = 12),
        strip.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        axis.text.x = element_text(size = 11),
        axis.text.y = element_text(size = 11),
        panel.grid.minor = element_blank(),
        panel.background = element_blank())

ggsave("figures/appendix/response_latency.pdf", width = 14, height = 6)

#How many respondents took the survey in English?----

#plot proportion of respondents who took survey in English per country

#summarize which language was used per country
lan_sum <- 
  all_countries_clean %>%
  group_by(country,UserLanguage) %>% 
  summarise(.,
            count = n()) %>%
  mutate(.,Language=UserLanguage)

#plot prop
ggplot(lan_sum, aes(x = country, y = count, fill = Language)) +
  geom_col(position = "fill") +
  xlab("Sample")+
  ylab("Proportion")
ggsave("figures/appendix/language_prop.pdf", width = 14, height = 6)


#Main analysis----

##Democratic Peace----

##Estimate Country Specific ATEs:

DP_bycountry <-
  all_countries_clean %>% 
  group_by(country) %>%
  do(tidy(lm_robust(outcome1_dem_new ~ 
                      DP_treat, data = .)))%>%
  filter(term!="(Intercept)")%>%
  mutate(.,
         meta = "Country Specific Replication",
         study = "Democratic Peace")%>%ungroup()%>%
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))



##Estimate Meta-Analytic Random Effects Model:
#DP_meta <- 
 # rma(yi=estimate, sei=std.error, data= DP_bycountry,
  #    slab = country)
#DP_meta$I2 #83.4

DP_meta <- 
  rma(yi=estimate, sei=std.error, data= DP_bycountry,
      slab = country) %>% broom::tidy(.,  conf.int=T) %>% 
  mutate(.,
         country = "All Countries", 
         meta = "Meta",
         study = "Democratic Peace",
         df = sum(DP_bycountry$df))


##Estimate original ATE from TW original data:
dem_peace_orig <- read_stata("data/original_experiments/TW_data_APSR.dta")
# Standardize outcome
dem_peace_orig <- 
  dem_peace_orig %>% 
  mutate(.,tw_support_attack = case_when( #scale from 1-5 instead of 0-4
           strike=="0"~1,
           strike=="1"~2,
           strike=="2"~3,
           strike=="3"~4,
           strike=="4"~5))

# Estimate model 
DP_orig <-
  dem_peace_orig %>% 
  lm_robust(tw_support_attack ~ 
              democ, data = .) %>% 
  tidy(.) %>% 
  filter(term == "democ") %>% 
  mutate(.,
         meta = "Original",
         country = "USA",
         study = "Democratic Peace")


##Audience Costs----

##Estimate Country Specific ATEs:

AC_bycountry <-
  all_countries_clean %>% 
  group_by(country) %>%
  do(tidy(lm_robust(outcome_aud_new ~ 
                      AC_treat, data = .)))%>%
  filter(term!="(Intercept)")%>%
  mutate(.,
         meta = "Country Specific Replication",
         study = "Audience Costs")%>%ungroup()%>%
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))


##Estimate Meta-Analytic Random Effects Model:
#AC_meta <- 
 # rma(yi=estimate, sei=std.error, data= AC_bycountry,
    #  slab = country)
#AC_meta$I2 #85.3

AC_meta <- 
  rma(yi=estimate, sei=std.error, data= AC_bycountry,
      slab = country) %>% tidy(.,  conf.int=T) %>% 
  mutate(.,
         country = "All Countries", 
         meta = "Meta",
         study = "Audience Costs",
         df = sum(AC_bycountry$df))

#Estimate ATE from BK original study:
audience_orig <- 
  get(load("data/original_experiments/Decomposing Data SSI.RData"))

#`Strategy` is treatment where (3 not_engage), (2 engage), (1 stay_out)
#We want to power on the traditional AC -- not_engage vs. stay out

audience_orig <- 
  audience_orig %>%
  mutate(.,
         treatment = case_when(
           Strategy == 3 ~ 1, 
           Strategy == 1 ~ 0),
         totApprov1 = case_when(
           totApprov== "-3"~1,
           totApprov== "-2"~2,
           totApprov== "-1"~3,
           totApprov== "0"~4,
           totApprov== "1"~5,
           totApprov== "2"~6,
           totApprov== "3"~7))

# Estimate model 
AC_orig <-
  audience_orig %>% 
  lm_robust(totApprov1~
              treatment,
            data = .) %>% 
  tidy(.) %>% 
  filter(term == "treatment") %>% 
  mutate(.,
         meta = "Original",
         country = "USA",
         study = "Audience Costs")

##International law----

##Estimate Country Specific ATEs:

IL_bycountry <-
  all_countries_clean %>% 
  group_by(country) %>%
  do(tidy(lm_robust(outcome_il_new ~
                      IL_treat, data = .)))%>%
  filter(term!="(Intercept)")%>%
  mutate(.,
         meta = "Country Specific Replication",
         study = "International Law")%>%ungroup()%>%
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))


##Estimate Meta-Analytic Random Effects Model:
#IL_meta <- 
 # rma(yi=estimate, sei=std.error, data= IL_bycountry,
  #    slab = country) 
#IL_meta$I2 #87.31

IL_meta <- 
  rma(yi=estimate, sei=std.error, data= IL_bycountry,
      slab = country) %>% tidy(.,  conf.int=T) %>% 
  mutate(.,
         country = "All Countries", 
         meta = "Meta",
         study = "International Law",
         df = sum(IL_bycountry$df))

##Estimate original ATE from Wallace data:
il_orig <- 
  read_dta('data/original_experiments/torture_law.dta')

# Estimate model
IL_orig <-
  il_orig %>% 
  lm_robust(torture7~
              treaty,
            data = .) %>% 
  tidy(.) %>% 
  filter(term == "treaty") %>% 
  mutate(.,
         meta = "Original",
         country = "USA",
         study = "International Law")


#IPE experiment----

##Estimate Country Specific ATEs:

IPE_bycountry <-
  all_countries_clean %>% 
  group_by(country) %>%
  do(tidy(lm_robust(outcome_ipe_new ~ 
                      IPE_treat, data = .)))%>%
  filter(term!="(Intercept)")%>%
  mutate(.,
         meta = "Country Specific Replication",
         study = "Reciprocity FDI")%>%ungroup()%>%
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))


##Estimate Meta-Analytic Random Effects Model:
#IPE_meta <- 
 # rma(yi=estimate, sei=std.error, data= IPE_bycountry,
  #    slab = country) 
#IPE_meta$I2 #97.6

IPE_meta <- 
  rma(yi=estimate, sei=std.error, data= IPE_bycountry,
      slab = country) %>% tidy(.,  conf.int=T) %>% 
  mutate(.,
         country = "All Countries", 
         meta = "Meta",
         study = "Reciprocity FDI",
         df = sum(IPE_bycountry$df))

##Estimate ATE from original IPE experiment data:
reciprocity_orig <-
  read_dta('data/original_experiments/reciprocityfollowup.dta')

# Clear status quo condition that we are not replicating 
reciprocity_orig <- 
  reciprocity_orig %>% 
  filter(.,
         current !=3) %>%
  mutate(.,
         outcome = case_when(
           pref == 0 ~ 1,
           pref == 0.25 ~ 2,
           pref == 0.50 ~ 3,
           pref == 0.75 ~ 4,
           pref == 1.00 ~ 5))

IPE_orig <-
  reciprocity_orig %>% 
  lm_robust(outcome~
              current,
            data = .) %>% 
  tidy(.) %>% 
  filter(term == "current") %>% 
  mutate(.,
         meta = "Original",
         country = "USA",
         study = "Reciprocity FDI")


#Combine different experiment data frames--------
all_experiments_df <-
  bind_rows(DP_bycountry,DP_meta,DP_orig,
            AC_bycountry,AC_meta,AC_orig,
            IL_bycountry,IL_meta,IL_orig,
            IPE_bycountry,IPE_meta,IPE_orig)%>%
  mutate(.,
         p_report = ifelse(meta=="Country Specific Replication",adjusted,p.value))

all_experiments_df$country <- 
  factor(all_experiments_df$country,
         levels = c("USA","Nigeria","Japan", "Israel", "India"
                    ,"Germany", "Brazil", "All Countries"))

all_experiments_df$meta <- 
  factor(all_experiments_df$meta,
         levels = c("Meta","Country Specific Replication","Original"))


ggplot(all_experiments_df, aes(x = as.numeric(estimate), y =as.factor(country), 
                               color =meta, shape = meta)) +
  geom_vline(xintercept = 0, color = "dodgerblue4", linetype = 2, size = 0.2) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high)) +
  geom_text(aes(label= paste(format(round(estimate, digits = 2), nsmall = 2))), hjust=0.3, vjust=-2, size = 3.3,
                  family = "Times") +
  geom_text(aes(label= paste("SE=",format(round(std.error, digits = 2), nsmall =2))), hjust=.3, vjust=-.7, size = 3.3, 
                     family = "Times")+
  #geom_text(aes(label= paste("p =", format(round(p_report, digits = 2), nsmall = 2))), hjust=0.3, vjust=-1.4, size = 3.3) +
  #geom_text(aes(label= paste("p. adj = ",format(round(adjusted, digits = 2), nsmall =2))), hjust=-.33, vjust=-0.5, size = 3.3, 
  #         family = "Times")+
  labs(x = "Main Effect",
       y = "") +
  facet_grid2(cols = vars(study), rows = vars(meta), scales= "free", space = "free_y")+
  scale_color_manual(values = c("grey40", "firebrick4", "dodgerblue4"))+
 # xlim(-2, 1.5)+
  theme_bw()+
  theme(strip.background =element_rect(fill="dodgerblue4"),
        strip.text = element_text(colour = 'white'),
        panel.grid.major = element_blank(),
        legend.position = "none",
        strip.text.x = element_text(size = 12),
        strip.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        axis.text.x = element_text(size = 11),
        axis.text.y = element_text(size = 11),
        panel.grid.minor = element_blank(),
        panel.background = element_blank())
ggsave("figures/main_text/main_meta_analysis_new.pdf", width = 13, height = 8.5)

#create the same figure but add demographic controls for appendix

##Democratic Peace----

##Estimate Country Specific ATEs:
DP_bycountry_cont <-
  all_countries_clean %>% 
  group_by(country) %>%
  do(tidy(lm_robust(outcome1_dem_new ~ 
      DP_treat + gender_fe + age + dem_norm + hawk + legal + 
      voting_fe + as.numeric(ideology) + university ##set covariates
    , data = .)))%>%
  filter(term!="(Intercept)")%>%
  mutate(.,
         meta = "Country Replication (w. Controls)",
         study = "Democratic Peace")%>%ungroup()%>%
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))

#keep only DP_treat estimates
DP_bycountry_cont2<-DP_bycountry_cont%>%
  filter(., term=="DP_treat")

##Estimate Meta-Analytic Random Effects Model:
DP_meta_cont <- 
  rma(yi=estimate, sei=std.error, data= DP_bycountry_cont2,
      slab = country) %>% broom::tidy(.,  conf.int=T) %>% 
  mutate(.,
         country = "All Countries", 
         meta = "Meta \n (w. Controls)",
         study = "Democratic Peace",
         df = sum(DP_bycountry_cont2$df))

##Audience Costs----

##Estimate Country Specific ATEs:
AC_bycountry_cont <-
  all_countries_clean %>% 
  group_by(country) %>%
  do(tidy(lm_robust(outcome_aud_new ~ 
      AC_treat  + gender_fe + age + dem_norm + hawk + legal + 
      voting_fe + as.numeric(ideology) + university ##set covariates
    , data = .)))%>%
  filter(term!="(Intercept)")%>%
  mutate(.,
         meta = "Country Replication (w. Controls)",
         study = "Audience Costs")%>%ungroup()%>%
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))

#keep only AC_treat estimates
AC_bycountry_cont2<-AC_bycountry_cont%>%
  filter(., term=="AC_treat")

##Estimate Meta-Analytic Random Effects Model:
AC_meta_cont <- 
  rma(yi=estimate, sei=std.error, data= AC_bycountry_cont2,
      slab = country) %>% tidy(.,  conf.int=T) %>% 
  mutate(.,
         country = "All Countries", 
         meta = "Meta \n (w. Controls)",
         study = "Audience Costs",
         df = sum(AC_bycountry_cont2$df))

##International law----

##Estimate Country Specific ATEs:

IL_bycountry_cont <-
  all_countries_clean %>% 
  group_by(country) %>%
  do(tidy(lm_robust(outcome_il_new ~
      IL_treat  + gender_fe + age + dem_norm + hawk + legal + 
      voting_fe + as.numeric(ideology) + university ##set covariates
    , data = .)))%>%
  filter(term!="(Intercept)")%>%
  mutate(.,
         meta = "Country Replication (w. Controls)",
         study = "International Law")%>%ungroup()%>%
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))

#keep only IL_treat estimates
IL_bycountry_cont2<-IL_bycountry_cont%>%
  filter(., term=="IL_treat")

##Estimate Meta-Analytic Random Effects Model:
IL_meta_cont <- 
  rma(yi=estimate, sei=std.error, data= IL_bycountry_cont2,
      slab = country) %>% tidy(.,  conf.int=T) %>% 
  mutate(.,
         country = "All Countries", 
         meta = "Meta \n (w. Controls)",
         study = "International Law",
         df = sum(IL_bycountry_cont2$df))

#IPE experiment----

##Estimate Country Specific ATEs:

IPE_bycountry_cont <-
  all_countries_clean %>% 
  group_by(country) %>%
  do(tidy(lm_robust(outcome_ipe_new ~ 
      IPE_treat + gender_fe + age + dem_norm + hawk + legal + 
      voting_fe + as.numeric(ideology) + university ##set covariates
    , data = .)))%>%
  filter(term!="(Intercept)")%>%
  mutate(.,
         meta = "Country Replication (w. Controls)",
         study = "Reciprocity FDI")%>%ungroup()%>%
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))

#keep only IPE_treat estimates
IPE_bycountry_cont2<-IPE_bycountry_cont%>%
  filter(., term=="IPE_treat")

##Estimate Meta-Analytic Random Effects Model:
IPE_meta_cont <- 
  rma(yi=estimate, sei=std.error, data= IPE_bycountry_cont2,
      slab = country) %>% tidy(.,  conf.int=T) %>% 
  mutate(.,
         country = "All Countries", 
         meta = "Meta \n (w. Controls)",
         study = "Reciprocity FDI",
         df = sum(IPE_bycountry_cont2$df))

#Combine different experiment data frames--------
all_experimentsCONT_df <-
  bind_rows(DP_bycountry_cont2,DP_meta_cont,
            AC_bycountry_cont2,AC_meta_cont,
            IL_bycountry_cont2,IL_meta_cont,
            IPE_bycountry_cont2,IPE_meta_cont)%>%
  mutate(.,
         p_report = ifelse(meta=="Country Replication (w. Controls)",adjusted,p.value))

all_experimentsCONT_df$country <- 
  factor(all_experimentsCONT_df$country,
         levels = c("USA","Nigeria","Japan", "Israel", "India"
                    ,"Germany", "Brazil", "All Countries"))

all_experimentsCONT_df$meta <- 
  factor(all_experimentsCONT_df$meta,
         levels = c("Meta \n (w. Controls)","Country Replication (w. Controls)"))


ggplot(all_experimentsCONT_df, aes(x = as.numeric(estimate), y =as.factor(country), 
                               color =meta, shape = meta)) +
  geom_vline(xintercept = 0, color = "dodgerblue4", linetype = 2, size = 0.2) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high)) +
  geom_text(aes(label= paste(format(round(estimate, digits = 2), nsmall = 2))), hjust=0.3, vjust=-2, size = 3.3,
            family = "Times") +
  geom_text(aes(label= paste("SE=",format(round(std.error, digits = 2), nsmall =2))), hjust=.3, vjust=-.7, size = 3.3, 
            family = "Times")+
  labs(x = "Main Effect",
       y = "") +
  facet_grid(cols = vars(study), rows = vars(meta), scales="free_y", space = "free")+
  scale_color_manual(values = c("grey40", "firebrick4", "dodgerblue4"))+
  xlim(-2, 1.5)+
  theme_bw()+
  theme(strip.background =element_rect(fill="dodgerblue4"),
        strip.text = element_text(colour = 'white'),
        panel.grid.major = element_blank(),
        legend.position = "none",
        strip.text.x = element_text(size = 12),
        strip.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        axis.text.x = element_text(size = 11),
        axis.text.y = element_text(size = 11),
        panel.grid.minor = element_blank(),
        panel.background = element_blank())
ggsave("figures/appendix/main_meta_analysis_cont.pdf", width = 8, height = 8)


#create appendix-table that corresponds to Figure 4

all_experiments_df$country_n<-ifelse(all_experiments_df$meta=="Original",
                                   "Original (USA)",as.character(all_experiments_df$country))

SE_dp<-subset(all_experiments_df$std.error,all_experiments_df$study=="Democratic Peace")
SE_ac<-subset(all_experiments_df$std.error,all_experiments_df$study=="Audience Costs")
SE_il<-subset(all_experiments_df$std.error,all_experiments_df$study=="International Law")
SE_ipe<-subset(all_experiments_df$std.error,all_experiments_df$study=="Reciprocity FDI")

Estimate_dp<-subset(all_experiments_df$estimate,all_experiments_df$study=="Democratic Peace")
Estimate_ac<-subset(all_experiments_df$estimate,all_experiments_df$study=="Audience Costs")
Estimate_il<-subset(all_experiments_df$estimate,all_experiments_df$study=="International Law")
Estimate_ipe<-subset(all_experiments_df$estimate,all_experiments_df$study=="Reciprocity FDI")

P_dp<-subset(all_experiments_df$p_report,all_experiments_df$study=="Democratic Peace")
P_ac<-subset(all_experiments_df$p_report,all_experiments_df$study=="Audience Costs")
P_il<-subset(all_experiments_df$p_report,all_experiments_df$study=="International Law")
P_ipe<-subset(all_experiments_df$p_report,all_experiments_df$study=="Reciprocity FDI")

N_dp<-subset(all_experiments_df$df,all_experiments_df$study=="Democratic Peace")
N_ac<-subset(all_experiments_df$df,all_experiments_df$study=="Audience Costs")
N_il<-subset(all_experiments_df$df,all_experiments_df$study=="International Law")
N_ipe<-subset(all_experiments_df$df,all_experiments_df$study=="Reciprocity FDI")

all_studies_wide<-data.frame(Samples = unique(all_experiments_df$country_n),
                             Estimate_dp = Estimate_dp,SE_dp = SE_dp,P_dp = P_dp,N_dp = N_dp,
                             Estimate_ac=Estimate_ac,SE_ac=SE_ac,P_ac=P_ac,N_ac=N_ac,
                             Estimate_il=Estimate_il, SE_il=SE_il,P_il=P_il,N_il=N_il,
                             Estimate_ipe=Estimate_ipe,SE_ipe=SE_ipe, P_ipe=P_ipe,N_ipe=N_ipe)

all_studies_wide<-rename(all_studies_wide,
                        "Estimate (DP)" = "Estimate_dp",
                        "Estimate (AC)" = "Estimate_ac",
                        "Estimate (IL)" = "Estimate_il",
                        "Estimate (FDI)" = "Estimate_ipe",
                        
                        "SE (DP)" = "SE_dp",
                        "SE (AC)" = "SE_ac",
                        "SE (IL)" = "SE_il",
                        "SE (FDI)" = "SE_ipe",
                        
                        "P value (DP)" = "P_dp",
                        "P value (AC)" = "P_ac",
                        "P value (IL)" = "P_il",
                        "P value (FDI)" = "P_ipe",
                        
                        "N (DP)" = "N_dp",
                        "N (AC)" = "N_ac",
                        "N (IL)" = "N_il",
                        "N (FDI)" = "N_ipe")


library(kableExtra)

write(kable(all_studies_wide, format = "latex", booktabs = T, caption = 
        "Meta analysis (Figure 3) in table form. \\label{tab:meta}", digits=2) %>% 
  add_header_above(c(" ", "Democratic Peace" = 4, "Audience Costs" = 4,
                     "International Law"=4, "Reciprocity (FDI)"=4)) %>% 
                     kable_styling(latex_options="scale_down"),'tables/meta_tab.tex')


#create appendix-table for meta figure with demographic controls:

all_experimentsCONT_df$country_n<-ifelse(all_experimentsCONT_df$meta=="Original",
                                     "Original (USA)",as.character(all_experimentsCONT_df$country))

SE_dp2<-subset(all_experimentsCONT_df$std.error,all_experimentsCONT_df$study=="Democratic Peace")
SE_ac2<-subset(all_experimentsCONT_df$std.error,all_experimentsCONT_df$study=="Audience Costs")
SE_il2<-subset(all_experimentsCONT_df$std.error,all_experimentsCONT_df$study=="International Law")
SE_ipe2<-subset(all_experimentsCONT_df$std.error,all_experimentsCONT_df$study=="Reciprocity FDI")

Estimate_dp2<-subset(all_experimentsCONT_df$estimate,all_experimentsCONT_df$study=="Democratic Peace")
Estimate_ac2<-subset(all_experimentsCONT_df$estimate,all_experimentsCONT_df$study=="Audience Costs")
Estimate_il2<-subset(all_experimentsCONT_df$estimate,all_experimentsCONT_df$study=="International Law")
Estimate_ipe2<-subset(all_experimentsCONT_df$estimate,all_experimentsCONT_df$study=="Reciprocity FDI")

P_dp2<-subset(all_experimentsCONT_df$p_report,all_experimentsCONT_df$study=="Democratic Peace")
P_ac2<-subset(all_experimentsCONT_df$p_report,all_experimentsCONT_df$study=="Audience Costs")
P_il2<-subset(all_experimentsCONT_df$p_report,all_experimentsCONT_df$study=="International Law")
P_ipe2<-subset(all_experimentsCONT_df$p_report,all_experimentsCONT_df$study=="Reciprocity FDI")

N_dp2<-subset(all_experimentsCONT_df$df,all_experimentsCONT_df$study=="Democratic Peace")
N_ac2<-subset(all_experimentsCONT_df$df,all_experimentsCONT_df$study=="Audience Costs")
N_il2<-subset(all_experimentsCONT_df$df,all_experimentsCONT_df$study=="International Law")
N_ipe2<-subset(all_experimentsCONT_df$df,all_experimentsCONT_df$study=="Reciprocity FDI")

all_studies_wide2<-data.frame(Samples = unique(all_experimentsCONT_df$country_n),
                             Estimate_dp2 = Estimate_dp2,SE_dp2 = SE_dp2,P_dp2 = P_dp2,N_dp2 = N_dp2,
                             Estimate_ac2=Estimate_ac2,SE_ac2=SE_ac2,P_ac2=P_ac2,N_ac2=N_ac2,
                             Estimate_il2=Estimate_il2, SE_il2=SE_il2,P_il2=P_il2,N_il2=N_il2,
                             Estimate_ipe2=Estimate_ipe2,SE_ipe2=SE_ipe2, P_ipe2=P_ipe2,N_ipe2=N_ipe2)

all_studies_wide2<-rename(all_studies_wide2,
                         "Estimate (DP)" = "Estimate_dp2",
                         "Estimate (AC)" = "Estimate_ac2",
                         "Estimate (IL)" = "Estimate_il2",
                         "Estimate (FDI)" = "Estimate_ipe2",
                         
                         "SE (DP)" = "SE_dp2",
                         "SE (AC)" = "SE_ac2",
                         "SE (IL)" = "SE_il2",
                         "SE (FDI)" = "SE_ipe2",
                         
                         "P value (DP)" = "P_dp2",
                         "P value (AC)" = "P_ac2",
                         "P value (IL)" = "P_il2",
                         "P value (FDI)" = "P_ipe2",
                         
                         "N (DP)" = "N_dp2",
                         "N (AC)" = "N_ac2",
                         "N (IL)" = "N_il2",
                         "N (FDI)" = "N_ipe2")

write(kable(all_studies_wide2, format = "latex", booktabs = T, caption = 
              "Meta analysis with controls (Figure A11) in table form. \\label{tab:meta2}", digits=2) %>% 
        add_header_above(c(" ", "Democratic Peace" = 4, "Audience Costs" = 4,
                           "International Law"=4, "Reciprocity (FDI)"=4)) %>% 
        kable_styling(latex_options="scale_down"),'tables/meta_tab2.tex')


# Implement Sign Generalization Test -----
us <- all_countries_clean %>% filter(country == "USA")
germany <- all_countries_clean %>% filter(country == "Germany")
india <- all_countries_clean %>% filter(country == "India")
israel <- all_countries_clean %>% filter(country == "Israel")
japan <- all_countries_clean %>% filter(country == "Japan")
nigeria <- all_countries_clean %>% filter(country == "Nigeria")
brazil <- all_countries_clean %>% filter(country == "Brazil")

country_list <- 
  list(us, germany, india, 
       israel, japan, nigeria, 
       brazil)

# Collect p values for each country testing whether we can reject the null of a positive effect (i.e. whether T-Pate is negative)
dem_sign_gen <- c()
for(i in 1:7){
  data_variation <- country_list[[i]]
  fit_variation <- lm_robust(outcome1_dem_new ~ DP_treat, data = data_variation)  # compute the SATE for each variation
  z <- summary(fit_variation)$coef["DP_treat", "Estimate"]/summary(fit_variation)$coef["DP_treat", "Std. Error"]
  dem_sign_gen[i] <- pnorm(z)  # one-sided p-value for each country's negative effect of democracy on support for conflict
}

# Implement Partial Conjunction Test
dp_sign_df <- 
  pct(dem_sign_gen) %>% 
  mutate(.,
         country = case_when(
           h_num == 1 ~ "us",
           h_num == 2 ~ "de",
           h_num == 3 ~ "in",
           h_num == 4 ~ "il",
           h_num == 5 ~ "jp",
           h_num == 6 ~ "ng",
           h_num == 7 ~ "br"), 
         study = "Democratic Peace")


## Sign Generalization Test AC------

## List of country dataset already exits "country_list"

# Collect p values for each country testing whether we can reject the null of a positive effect (i.e. whether T-Pate is negative)
ac_sign_gen <- c()
for(i in 1:7){
  data_variation <- country_list[[i]]
  fit_variation <- lm_robust(outcome_aud_new ~ AC_treat, data = data_variation)  # compute the SATE for each variation
  z <- summary(fit_variation)$coef["AC_treat", "Estimate"]/summary(fit_variation)$coef["AC_treat", "Std. Error"]
  ac_sign_gen[i] <- pnorm(z)  # one-sided p-value for each country's negative effect of democracy on support for conflict
}

# Implement Partial Conjunction Test
ac_sign_df <- 
  pct(ac_sign_gen) %>% 
  mutate(.,
         country = case_when(
           h_num == 1 ~ "us",
           h_num == 2 ~ "de",
           h_num == 3 ~ "in",
           h_num == 4 ~ "il",
           h_num == 5 ~ "jp",
           h_num == 6 ~ "ng",
           h_num == 7 ~ "br"), 
         study = "Audience Cost")

## Sign Generalization Test IL------

## List of country dataset already exits "country_list"

# Collect p values for each country testing whether we can reject the null of a positive effect (i.e. whether T-Pate is negative)
il_sign_gen <- c()
for(i in 1:7){
  data_variation <- country_list[[i]]
  fit_variation <- lm_robust(outcome_il_new ~ IL_treat, data = data_variation)  # compute the SATE for each variation
  z <- summary(fit_variation)$coef["IL_treat", "Estimate"]/summary(fit_variation)$coef["IL_treat", "Std. Error"]
  il_sign_gen[i] <- pnorm(z)  # one-sided p-value for each country's negative effect of democracy on support for conflict
}

# Implement Partial Conjunction Test
il_sign_df <- 
  pct(il_sign_gen) %>% 
  mutate(.,
         country = case_when(
           h_num == 1 ~ "us",
           h_num == 2 ~ "de",
           h_num == 3 ~ "in",
           h_num == 4 ~ "il",
           h_num == 5 ~ "jp",
           h_num == 6 ~ "ng",
           h_num == 7 ~ "br"), 
         study = "International Law")


# Sign Generalization Test IPE------

## List of country dataset already exits "country_list"

# Collect p values for each country testing whether we can reject the null of a !!!!negative!!!! effect (i.e. whether T-Pate is Positive)
ipe_sign_gen <- c()
for(i in 1:7){
  data_variation <- country_list[[i]]
  fit_variation <- lm_robust(outcome_ipe_new ~ IPE_treat, data = data_variation)  # compute the SATE for each variation
  z <- summary(fit_variation)$coef["IPE_treat", "Estimate"]/summary(fit_variation)$coef["IPE_treat", "Std. Error"]
  ipe_sign_gen[i] <- 1- pnorm(z)  # one-sided p-value for each country's negative effect of democracy on support for conflict
}

# Implement Partial Conjunction Test
ipe_sign_df <- 
  pct(ipe_sign_gen) %>% 
  mutate(.,
         country = case_when(
           h_num == 1 ~ "us",
           h_num == 2 ~ "de",
           h_num == 3 ~ "in",
           h_num == 4 ~ "il",
           h_num == 5 ~ "jp",
           h_num == 6 ~ "ng",
           h_num == 7 ~ "br"), 
         study = "Reciprocity FDI")


###Combine all sign generalization tests----
sign_general_combine <-
  bind_rows(dp_sign_df, ac_sign_df, 
            il_sign_df, ipe_sign_df)

# Plot P values from Partial Conjunction Test
ggplot(sign_general_combine, aes(threshold, p_value, country = country))+
  geom_flag()+
  ylim(0,.75)+# Maybe change Ylim with real data where p values will be much smaller
  labs(y="Partial Conjunction P-Value",
       x = "Threshold")+
  geom_text(aes(label= sprintf("%.2f", round(p_value,digits =2)), 
              #  position = position_dodge(width = .6),
                vjust = -2.8,
              #  size = 3
              ))+
  scale_x_continuous(breaks = c(1:7,1))+
  geom_hline(yintercept=0.05, linetype="dashed", color = "dodgerblue4")+
  facet_wrap(~study)+
  #xlim(-0.05, 0.05)+
  theme_bw()+
  theme(strip.background =element_rect(fill="dodgerblue4"),
        strip.text = element_text(colour = 'white'),
        panel.grid.major = element_blank(),
        legend.position = "none",
        strip.text.x = element_text(size = 12),
        strip.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        axis.text.x = element_text(size = 11),
        axis.text.y = element_text(size = 11),
        panel.grid.minor = element_blank(),
        panel.background = element_blank())
ggsave("figures/main_text/sign_generalization_main.pdf", width = 8, height = 7)

#create appendix-table of the sign generalization figure

sign_general_combine_wide <- 
  sign_general_combine%>%
  mutate(.,
         country2 = case_when(
           country == "us"~"USA",
           country == "de"~"Germany",
           country == "in"~"India",
           country == "il"~"Israel",
           country == "jp"~"Japan",
           country == "ng"~"Nigeria",
           country == "br"~"Brazil")) 

sign_general_combine_wide<- sign_general_combine_wide[order(sign_general_combine_wide$study,
                                                            sign_general_combine_wide$country2),] 

thresh_dp<-subset(sign_general_combine_wide$threshold,sign_general_combine_wide$study=="Democratic Peace")
thresh_ac<-subset(sign_general_combine_wide$threshold,sign_general_combine_wide$study=="Audience Cost")
thresh_il<-subset(sign_general_combine_wide$threshold,sign_general_combine_wide$study=="International Law")
thresh_ipe<-subset(sign_general_combine_wide$threshold,sign_general_combine_wide$study=="Reciprocity FDI")

p_dp<-subset(sign_general_combine_wide$p_value,sign_general_combine_wide$study=="Democratic Peace")
p_ac<-subset(sign_general_combine_wide$p_value,sign_general_combine_wide$study=="Audience Cost")
p_il<-subset(sign_general_combine_wide$p_value,sign_general_combine_wide$study=="International Law")
p_ipe<-subset(sign_general_combine_wide$p_value,sign_general_combine_wide$study=="Reciprocity FDI")

sign_general_combine_wide<-data.frame(Samples = unique(sign_general_combine_wide$country2),
                             thresh_dp = thresh_dp,p_dp = p_dp,
                             thresh_ac=thresh_ac,p_ac=p_ac,
                             thresh_il=thresh_il, p_il=p_il,
                             thresh_ipe=thresh_ipe,p_ipe=p_ipe)

sign_general_combine_wide<-rename(sign_general_combine_wide,
                         "Threshold (DP)" = "thresh_dp",
                         "Threshold (AC)" = "thresh_ac",
                         "Threshold (IL)" = "thresh_il",
                         "Threshold (FDI)" = "thresh_ipe",
                         
                         "P Value (DP)" = "p_dp",
                         "P Value (AC)" = "p_ac",
                         "P Value (IL)" = "p_il",
                         "P Value (FDI)" = "p_ipe")

write(kable(sign_general_combine_wide, format = "latex", booktabs = T, caption = 
              "Sign Generalization (Figure 4) in table form. \\label{tab:sign}", digits=2) %>% 
        add_header_above(c(" ", "Democratic Peace" = 2, "Audience Costs" = 2,
                           "International Law"=2, "Reciprocity (FDI)"=2)) %>% 
        kable_styling(latex_options="scale_down"),'tables/sign_tab.tex')


#Plausibility by experiment----

#DP
plaus_DP<-all_countries_clean%>%
  dplyr::select(country,plausible_DP1)%>%
  mutate(.,
         experiment = "Democratic Peace")%>%
  rename(plausible = plausible_DP1)%>% drop_na() 
#AC
plaus_AC<-all_countries_clean%>%
  dplyr::select(country,plausible_AC1)%>%
  mutate(.,
         experiment = "Audience Costs")%>%
  rename(plausible = plausible_AC1)%>% drop_na()
#IL
plaus_IL<-all_countries_clean%>%
  dplyr::select(country,plausible_IL1)%>%
  mutate(.,
         experiment = "International Law")%>%
  rename(plausible = plausible_IL1)%>% drop_na()
#IPE
plaus_IPE<-all_countries_clean%>%
  dplyr::select(country,plausible_IPE1)%>%
  mutate(.,
         experiment = "Receprocity FDI")%>%
  rename(plausible = plausible_IPE1)%>% drop_na()

plaus_all<-rbind(plaus_DP,plaus_AC,
                 plaus_IL,plaus_IPE)


ggplot(plaus_all, aes(x= plausible,  group=country)) + 
  geom_bar(aes(y = ..prop..), stat="count") +
  geom_text(aes(label = scales::percent(..prop..),
                 y= ..prop.. ), stat= "count", vjust = -.5, size = 2) +
  labs(y = "Percent", x="Scenario plausible") +
  facet_grid(country~experiment) +
  theme_bw()+
  theme(text = element_text(size = 14, family = "Times"),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank()
  )+
  scale_y_continuous(labels = scales::percent, limits = c(0,1))

ggsave("figures/appendix/plaus.pdf",width = 14, height = 7)


#Manipulation checks: Estimate ATE on treatment recall-------
#DP
us_dp_mnp <- lm_robust(DP_man ~ DP_treat, subset = country == "USA", data = all_countries_clean)
br_dp_mnp <- lm_robust(DP_man ~ DP_treat, subset = country == "Brazil", data = all_countries_clean)
gr_dp_mnp <- lm_robust(DP_man ~ DP_treat, subset = country == "Germany", data = all_countries_clean)
in_dp_mnp <- lm_robust(DP_man ~ DP_treat, subset = country == "India", data = all_countries_clean)
il_dp_mnp <- lm_robust(DP_man ~ DP_treat, subset = country == "Israel", data = all_countries_clean)
jp_dp_mnp <- lm_robust(DP_man ~ DP_treat, subset = country == "Japan", data = all_countries_clean)
ng_dp_mnp <- lm_robust(DP_man ~ DP_treat, subset = country == "Nigeria", data = all_countries_clean)

#AC
us_ac_mnp <- lm_robust(AC_man ~ AC_treat, subset = country == "USA", data = all_countries_clean)
br_ac_mnp <- lm_robust(AC_man ~ AC_treat, subset = country == "Brazil", data = all_countries_clean)
gr_ac_mnp <- lm_robust(AC_man ~ AC_treat, subset = country == "Germany", data = all_countries_clean)
in_ac_mnp <- lm_robust(AC_man ~ AC_treat, subset = country == "India", data = all_countries_clean)
il_ac_mnp <- lm_robust(AC_man ~ AC_treat, subset = country == "Israel", data = all_countries_clean)
jp_ac_mnp <- lm_robust(AC_man ~ AC_treat, subset = country == "Japan", data = all_countries_clean)
ng_ac_mnp <- lm_robust(AC_man ~ AC_treat, subset = country == "Nigeria", data = all_countries_clean)


#IL
us_il_mnp <- lm_robust(IL_man ~ IL_treat, subset = country == "USA", data = all_countries_clean)
br_il_mnp <- lm_robust(IL_man ~ IL_treat, subset = country == "Brazil", data = all_countries_clean)
gr_il_mnp <- lm_robust(IL_man ~ IL_treat, subset = country == "Germany", data = all_countries_clean)
in_il_mnp <- lm_robust(IL_man ~ IL_treat, subset = country == "India", data = all_countries_clean)
il_il_mnp <- lm_robust(IL_man ~ IL_treat, subset = country == "Israel", data = all_countries_clean)
jp_il_mnp <- lm_robust(IL_man ~ IL_treat, subset = country == "Japan", data = all_countries_clean)
ng_il_mnp <- lm_robust(IL_man ~ IL_treat, subset = country == "Nigeria", data = all_countries_clean)


#IPE
us_ip_mnp <- lm_robust(IPE_man ~ IPE_treat, subset = country == "USA", data = all_countries_clean)
br_ip_mnp <- lm_robust(IPE_man ~ IPE_treat, subset = country == "Brazil", data = all_countries_clean)
gr_ip_mnp <- lm_robust(IPE_man ~ IPE_treat, subset = country == "Germany", data = all_countries_clean)
in_ip_mnp <- lm_robust(IPE_man ~ IPE_treat, subset = country == "India", data = all_countries_clean)
il_ip_mnp <- lm_robust(IPE_man ~ IPE_treat, subset = country == "Israel", data = all_countries_clean)
jp_ip_mnp <- lm_robust(IPE_man ~ IPE_treat, subset = country == "Japan", data = all_countries_clean)
ng_ip_mnp <- lm_robust(IPE_man ~ IPE_treat, subset = country == "Nigeria", data = all_countries_clean)

#print table
texreg(list(
  us_dp_mnp,br_dp_mnp,gr_dp_mnp,in_dp_mnp,il_dp_mnp,jp_dp_mnp,ng_dp_mnp,
  us_ac_mnp,br_ac_mnp,gr_ac_mnp,in_ac_mnp,il_ac_mnp,jp_ac_mnp,ng_ac_mnp,
  us_il_mnp,br_il_mnp,gr_il_mnp,in_il_mnp,il_il_mnp,jp_il_mnp,ng_il_mnp,
  us_ip_mnp,br_ip_mnp,gr_ip_mnp,in_ip_mnp,il_ip_mnp,jp_ip_mnp,ng_ip_mnp),
       label = "tab:treat_manip",
       caption.above = T,
       no.margin = T,
       include.ci = FALSE,
       custom.header = list("Country is Democracy" = 1:7,
                            "Leader back down" = 8:14,
                            "Torture Violates Law" = 15:21,
                            "Investment Made Harder" = 22:28),
       custom.model.names	= c("USA", "BRZ", "GRM", "IND", "ISL", "JPN", "NGR",
                              "USA", "BRZ", "GRM", "IND", "ISL", "JPN", "NGR",
                              "USA", "BRZ", "GRM", "IND", "ISL", "JPN", "NGR",
                              "USA", "BRZ", "GRM", "IND", "ISL", "JPN", "NGR"),
       custom.coef.map = list("DP_treat" = "Democracy", "AC_treat" = "Engage",  
                              "IL_treat" = "IL Law", "IPE_treat" = "Harder Invest"),
       fontsize = "tiny",
      stars = 0.05,
       digits = 2,
       include.rsquared = F, include.rmse = F,
       caption = "Manipulation Test: Treatment Effects on Correct Recall",
       font.size = "scriptsize",
       table = T,
       file = "tables/manipulation_check.tex")


#Heterogenous treatment effects ------

##Plot moderators descriptively-----
#dem_norm
mod_dem<-all_countries_clean%>%
  dplyr::select(country,dem_norm)%>%
  mutate(.,
         moderator = "Democratic Norms")%>%
  rename(value = dem_norm) %>% drop_na()

#hawk
mod_hawk<-all_countries_clean%>%
  dplyr::select(country,hawk)%>%
  mutate(.,
         moderator = "Hawkishness")%>%
  rename(value = hawk) %>% drop_na()

#legal
mod_legal<-all_countries_clean%>%
  dplyr::select(country,legal)%>%
  mutate(.,
         moderator = "Intl Legal Obligation")%>%
  rename(value = legal)%>% drop_na()

#bind
moderators<-rbind(mod_dem,mod_hawk,
                 mod_legal) 


summary.moderators <- moderators %>% 
  group_by(country,moderator) %>% 
  summarise(mean.value = mean(value))

#plot
ggplot(moderators, aes(x=value)) + 
  geom_histogram(fill = "dodgerblue4")+
  labs(x = "Moderator Distribution",
       y = "Count")+
  geom_vline(data = summary.moderators, aes(xintercept = mean.value), colour = "red") +
  facet_grid(country~moderator)+
  theme_bw()+
  theme(text = element_text(size = 14, family = "Times"),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank()
  )
ggsave("figures/appendix/moderators.pdf",width = 6, height = 6)


## Estimate HTE for DP------
all_countries_clean<-as.data.frame(all_countries_clean)

interflex(Y = "outcome1_dem_new", D = "DP_treat", X = "dem_norm",
          data = all_countries_clean,
          estimator = "raw", 
          na.rm = T,
          main = "Marginal Effects")


mod_norms<- interflex(Y = "outcome1_dem_new", D = "DP_treat", X = "dem_norm", 
                 Z = c("country_fe", "gender_fe", "age","education_fe", "voting_fe"), data = all_countries_clean,
                 estimator = "binning", vcov.type = "robust",
                 na.rm = T,
                 main = "Marginal Effects")


plot(mod_norms,
     theme.bw = T, show.grid = F,
     Ylabel = "Support for Attack", 
     Dlabel = "Democracy", Xlabel = "Support for Democratic Norms",
     cex.lab = .5)

ggsave("figures/appendix/hte_dp.pdf",width = 6, height = 4)


## Estimate HTE for AC------

interflex(Y = "outcome_aud_new", D = "AC_treat", X = "hawk", 
          data = all_countries_clean,
          estimator = "raw", 
          na.rm = T,
          main = "Marginal Effects")
mod_hawk<- interflex(Y = "outcome_aud_new", D = "AC_treat", X = "hawk", 
                      Z = c("country_fe", "gender_fe", "age","education_fe", "voting_fe"), data = all_countries_clean,
                      estimator = "binning", vcov.type = "robust",
                      na.rm = T,
                      main = "Marginal Effects")

plot(mod_hawk,
     theme.bw = T, show.grid = F,
     Ylabel = "Leader Approval", 
     Dlabel = "Back Down", Xlabel = "Hawkishness",
     cex.lab = .5)

ggsave("figures/appendix/hte_ac.pdf",width = 6, height = 4)


## Estimate HTE for IL------

interflex(Y = "outcome_il_new", D = "IL_treat", X = "legal", data = all_countries_clean,
          estimator = "raw", 
          na.rm = T,
          main = "Marginal Effects")

mod_leg<- interflex(Y = "outcome_il_new", D = "IL_treat", X = "legal", 
                     Z = c("country_fe", "gender_fe", "age","education_fe", "voting_fe"), data = all_countries_clean,
                     estimator = "binning", vcov.type = "robust",
                     na.rm = T,
                     main = "Marginal Effects")

plot(mod_leg,
     theme.bw = T, show.grid = F,
     Ylabel = "Support Torture", 
     Dlabel = "International Law", Xlabel = "Legal Obligation",
     cex.lab = .5)

ggsave("figures/appendix/hte_il.pdf",width = 8, height = 4)



##Create HTE Table------
dp_hte_tbl <- 
lm_lin(outcome1_dem_new~DP_treat,
       covariates = ~dem_norm+country_fe+gender_fe+age+education_fe+voting_fe,
       data = all_countries_clean) 

ac_hte_tbl <- 
  lm_lin(outcome_aud_new~AC_treat,
         covariates = ~hawk+country_fe+gender_fe+age+education_fe+voting_fe,
         data = all_countries_clean) 


il_hte_tbl <- 
  lm_lin(outcome_il_new~IL_treat,
         covariates = ~legal+country_fe+gender_fe+age+education_fe+voting_fe,
         data = all_countries_clean) 

texreg(list(
  dp_hte_tbl, ac_hte_tbl, il_hte_tbl),
  label = "tab:hte_table",
  caption.above = T,
  no.margin = T,
  include.ci = FALSE,
  custom.header = list("Democratic Peace" = 1,
                       "Audience Costs" = 2,
                       "International Law" = 3),
  custom.coef.map = list("DP_treat" = "Democracy", "dem_norm_c" = "Dem Norms",  
                         "DP_treat:dem_norm_c" = "Dem*Norms",
                         "AC_treat" = "Back Down", "hawk_c" = "Hawkish", 
                         "AC_treat:hawk_c" = "BD*Hawk",
                         "IL_treat" = "Intl Law", "legal_c" = "Legal Obligation", "IL_treat:legal_c" = "IL*Oblig"),
  stars = 0.05,
  digits = 3,
  include.rsquared = F, include.rmse = F,
  caption = "Moderation Tests",
  font.size = "tiny",
  table = T,
  custom.note = "%stars. Regressions interact treatment with covariates (gender, age, education, voting status, country).",
  file = "tables/hte_table.tex")


#create a HTE table for each country country
#USA
us_new<-all_countries_clean%>%
  filter(.,country=="USA")

dp_hte_USA <- 
  lm_lin(outcome1_dem_new~DP_treat,
         covariates = ~dem_norm+gender_fe+age+education_fe+voting_fe,
         data = us_new) 
summary(dp_hte_USA)

ac_hte_USA <- lm_lin(outcome_aud_new~AC_treat,
       covariates = ~hawk+gender_fe+age+education_fe+voting_fe,
       data = us_new) 

summary(ac_hte_USA)

il_hte_USA <- 
  lm_lin(outcome_il_new~IL_treat,
       covariates = ~legal+gender_fe+age+education_fe+voting_fe,
       data = us_new) 
summary(il_hte_USA)

#Nigeria
nig_new<-all_countries_clean%>%
  filter(.,country=="Nigeria")

dp_hte_nig <- 
  lm_lin(outcome1_dem_new~DP_treat,
         covariates = ~dem_norm+gender_fe+age+voting_fe,
         data = nig_new) 
summary(dp_hte_nig)

ac_hte_nig <- lm_lin(outcome_aud_new~AC_treat,
                     covariates = ~hawk+gender_fe+age+voting_fe,
                     data = nig_new) 

summary(ac_hte_nig)

il_hte_nig <- 
  lm_lin(outcome_il_new~IL_treat,
         covariates = ~legal+gender_fe+age+voting_fe,
         data = nig_new) 
summary(il_hte_nig)


#Japan
jap_new<-all_countries_clean%>%
  filter(.,country=="Japan")

dp_hte_jap <- 
  lm_lin(outcome1_dem_new~DP_treat,
         covariates = ~dem_norm+gender_fe+age+education_fe+voting_fe,
         data = jap_new) 
summary(dp_hte_jap)

ac_hte_jap <- lm_lin(outcome_aud_new~AC_treat,
                     covariates = ~hawk+gender_fe+age+education_fe+voting_fe,
                     data = jap_new) 

summary(ac_hte_jap)

il_hte_jap <- 
  lm_lin(outcome_il_new~IL_treat,
         covariates = ~legal+gender_fe+age+education_fe+voting_fe,
         data = jap_new) 
summary(il_hte_jap)


#India
india_new<-all_countries_clean%>%
  filter(.,country=="India")

dp_hte_ind <- 
  lm_lin(outcome1_dem_new~DP_treat,
         covariates = ~dem_norm+gender_fe+age+education_fe+voting_fe,
         data = india_new) 
summary(dp_hte_ind)

ac_hte_ind <- lm_lin(outcome_aud_new~AC_treat,
                     covariates = ~hawk+gender_fe+age+education_fe+voting_fe,
                     data = india_new) 

summary(ac_hte_ind)

il_hte_ind <- 
  lm_lin(outcome_il_new~IL_treat,
         covariates = ~legal+gender_fe+age+education_fe+voting_fe,
         data = india_new) 
summary(il_hte_ind)


#Israel
israel_new<-all_countries_clean%>%
  filter(.,country=="Israel")


dp_hte_isr <- 
  lm_lin(outcome1_dem_new~DP_treat,
         covariates = ~dem_norm+gender_fe+age+education_fe+voting_fe,
         data = israel_new) 
summary(dp_hte_isr)

ac_hte_isr <- lm_lin(outcome_aud_new~AC_treat,
                     covariates = ~hawk+gender_fe+age+education_fe+voting_fe,
                     data = israel_new) 

summary(ac_hte_isr)

il_hte_isr <- 
  lm_lin(outcome_il_new~IL_treat,
         covariates = ~legal+gender_fe+age+education_fe+voting_fe,
         data = israel_new) 
summary(il_hte_isr)


#Germany
germany_new<-all_countries_clean%>%
  filter(.,country=="Germany")


dp_hte_ger <- 
  lm_lin(outcome1_dem_new~DP_treat,
         covariates = ~dem_norm+gender_fe+age+education_fe+voting_fe,
         data = germany_new) 
summary(dp_hte_ger)

ac_hte_ger <- lm_lin(outcome_aud_new~AC_treat,
                     covariates = ~hawk+gender_fe+age+education_fe+voting_fe,
                     data = germany_new) 

summary(ac_hte_ger)

il_hte_ger <- 
  lm_lin(outcome_il_new~IL_treat,
         covariates = ~legal+gender_fe+age+education_fe+voting_fe,
         data = germany_new) 
summary(il_hte_ger)



#Brazil
brazil_new<-all_countries_clean%>%
  filter(.,country=="Brazil")


dp_hte_Br <- 
  lm_lin(outcome1_dem_new~DP_treat,
         covariates = ~dem_norm+gender_fe+age+education_fe+voting_fe,
         data = brazil_new) 
summary(dp_hte_Br)

ac_hte_Br <- lm_lin(outcome_aud_new~AC_treat,
                     covariates = ~hawk+gender_fe+age+education_fe+voting_fe,
                     data = brazil_new) 

summary(ac_hte_Br)

il_hte_Br <- 
  lm_lin(outcome_il_new~IL_treat,
         covariates = ~legal+gender_fe+age+education_fe+voting_fe,
         data = brazil_new) 
summary(il_hte_Br)



#DP By country
DP_hte<-texreg(list(
  dp_hte_Br, dp_hte_ger, dp_hte_ind,dp_hte_isr,dp_hte_jap,dp_hte_nig,dp_hte_USA),
  label = "tab:hte_DP",
  caption.above = T,
  no.margin = T,
  include.ci = FALSE,
  custom.model.names = c("Brazil", "Germany", "India","Israel","Japan","Nigeria","USA"),
  custom.coef.map = list("DP_treat" = "Democracy", "dem_norm_c" = "Dem Norms",  
                         "DP_treat:dem_norm_c" = "Dem*Norms"),
  stars = 0.05,
  digits = 3,
  include.rsquared = F, include.rmse = F,
  caption = "Moderation Test (Democratic Peace)",
  font.size = "tiny",
  table = T)

DP_hte <- gsub("\\begin{tabular}"
               ,"\\resizebox{0.8\\textwidth}{!}{\n\\begin{tabular}"
               ,DP_hte
               ,fixed=T)

DP_hte <- gsub("\\end{tabular}"
           ,"\\end{tabular}}"
           ,DP_hte
           ,fixed=T)

write_file(DP_hte, 'tables/hte_DP.tex')


#AC By country
AC_hte<-texreg(list(
  ac_hte_Br, ac_hte_ger, ac_hte_ind,ac_hte_isr,ac_hte_jap,ac_hte_nig,ac_hte_USA),
  label = "tab:hte_AC",
  caption.above = T,
  no.margin = T,
  include.ci = FALSE,
  custom.model.names = c("Brazil", "Germany", "India","Israel","Japan","Nigeria","USA"),
  custom.coef.map = list("AC_treat" = "Back Down", "hawk_c" = "Hawkish", 
                         "AC_treat:hawk_c" = "BD*Hawk"),
  stars = 0.05,
  digits = 3,
  include.rsquared = F, include.rmse = F,
  caption = "Moderation Test (Audience Costs)",
  font.size = "tiny",
  table = T)

AC_hte <- gsub("\\begin{tabular}"
               ,"\\resizebox{0.8\\textwidth}{!}{\n\\begin{tabular}"
               ,AC_hte
               ,fixed=T)

AC_hte <- gsub("\\end{tabular}"
               ,"\\end{tabular}}"
               ,AC_hte
               ,fixed=T)

write_file(AC_hte, 'tables/hte_AC.tex')


#IL By country
IL_hte<-texreg(list(
  il_hte_Br, il_hte_ger, il_hte_ind,il_hte_isr,il_hte_jap,il_hte_nig,il_hte_USA),
  label = "tab:hte_IL",
  caption.above = T,
  no.margin = T,
  include.ci = FALSE,
  custom.model.names = c("Brazil", "Germany", "India","Israel","Japan","Nigeria","USA"),
  custom.coef.map = list("IL_treat" = "Intl Law", "legal_c" = "Legal Obligation", 
                         "IL_treat:legal_c" = "IL*Oblig"),
  stars = 0.05,
  digits = 3,
  include.rsquared = F, include.rmse = F,
  caption = "Moderation Test (International Law)",
  font.size = "tiny",
  table = T)

IL_hte <- gsub("\\begin{tabular}"
               ,"\\resizebox{0.8\\textwidth}{!}{\n\\begin{tabular}"
               ,IL_hte
               ,fixed=T)

IL_hte <- gsub("\\end{tabular}"
               ,"\\end{tabular}}"
               ,IL_hte
               ,fixed=T)

write_file(IL_hte, 'tables/hte_IL.tex')


## Create HTE figures for each country-study pair ------

#USA
us_dp_hte<- interflex(Y = "outcome1_dem_new", D = "DP_treat", X = "dem_norm", 
                      Z = c("gender_fe", "age","education_fe", "voting_fe"), data = us_new,
                      estimator = "binning", vcov.type = "robust",
                      na.rm = T,
                      main = "Marginal Effects")

us_ac_hte<- interflex(Y = "outcome_aud_new", D = "AC_treat", X = "hawk", 
                     Z = c("gender_fe", "age","education_fe", "voting_fe"), data = us_new,
                     estimator = "binning", vcov.type = "robust",
                     na.rm = T,
                     main = "Marginal Effects")

us_il_hte<- interflex(Y = "outcome_il_new", D = "IL_treat", X = "legal", 
                    Z = c("gender_fe", "age","education_fe", "voting_fe"), data = us_new,
                    estimator = "binning", vcov.type = "robust",
                    na.rm = T,
                    main = "Marginal Effects")


#Nigeria
nig_dp_hte<- interflex(Y = "outcome1_dem_new", D = "DP_treat", X = "dem_norm", 
                      Z = c("gender_fe", "age","education_fe", "voting_fe"), data = nig_new,
                      estimator = "binning", vcov.type = "robust",
                      na.rm = T,
                      main = "Marginal Effects")

nig_ac_hte<- interflex(Y = "outcome_aud_new", D = "AC_treat", X = "hawk", 
                      Z = c("gender_fe", "age","education_fe", "voting_fe"), data = nig_new,
                      estimator = "binning", vcov.type = "robust",
                      na.rm = T,
                      main = "Marginal Effects")

nig_il_hte<- interflex(Y = "outcome_il_new", D = "IL_treat", X = "legal", 
                      Z = c("gender_fe", "age","education_fe", "voting_fe"), data = nig_new,
                      estimator = "binning", vcov.type = "robust",
                      na.rm = T,
                      main = "Marginal Effects")

#Japan
jap_dp_hte<- interflex(Y = "outcome1_dem_new", D = "DP_treat", X = "dem_norm", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = jap_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

jap_ac_hte<- interflex(Y = "outcome_aud_new", D = "AC_treat", X = "hawk", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = jap_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

jap_il_hte<- interflex(Y = "outcome_il_new", D = "IL_treat", X = "legal", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = jap_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

#Israel
isr_dp_hte<- interflex(Y = "outcome1_dem_new", D = "DP_treat", X = "dem_norm", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = israel_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

isr_ac_hte<- interflex(Y = "outcome_aud_new", D = "AC_treat", X = "hawk", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = israel_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

isr_il_hte<- interflex(Y = "outcome_il_new", D = "IL_treat", X = "legal", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = israel_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")
#India
ind_dp_hte<- interflex(Y = "outcome1_dem_new", D = "DP_treat", X = "dem_norm", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = india_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

ind_ac_hte<- interflex(Y = "outcome_aud_new", D = "AC_treat", X = "hawk", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = india_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

ind_il_hte<- interflex(Y = "outcome_il_new", D = "IL_treat", X = "legal", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = india_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

#Germany
ger_dp_hte<- interflex(Y = "outcome1_dem_new", D = "DP_treat", X = "dem_norm", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = germany_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

ger_ac_hte<- interflex(Y = "outcome_aud_new", D = "AC_treat", X = "hawk", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = germany_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

ger_il_hte<- interflex(Y = "outcome_il_new", D = "IL_treat", X = "legal", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = germany_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")
#Brazil
br_dp_hte<- interflex(Y = "outcome1_dem_new", D = "DP_treat", X = "dem_norm", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = brazil_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

br_ac_hte<- interflex(Y = "outcome_aud_new", D = "AC_treat", X = "hawk", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = brazil_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

br_il_hte<- interflex(Y = "outcome_il_new", D = "IL_treat", X = "legal", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = brazil_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

library(patchwork)

#plot

#Brazil
pbr_dp_hte<-plot(br_dp_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support for Attack", 
     Dlabel = "Democracy", Xlabel = "Support for Democratic Norms",
     cex.lab = .5) +
  ggtitle('Democratic Peace (Brazil)')

pbr_ac_hte<-plot(br_ac_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Leader Approval", 
     Dlabel = "Back Down", Xlabel = "Hawkishness",
     cex.lab = .5) +
  ggtitle("Audience Costs (Brazil)")

pbr_il_hte<-plot(br_il_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support Torture", 
     Dlabel = "International Law", Xlabel = "Legal Obligation",
     cex.lab = .5) +
  ggtitle("International Law (Brazil)")

#Germany
pger_dp_hte<-plot(ger_dp_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support for Attack", 
     Dlabel = "Democracy", Xlabel = "Support for Democratic Norms",
     cex.lab = .5)+
  ggtitle("Democratic Peace (Germany)")

pger_ac_hte<-plot(ger_ac_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Leader Approval", 
     Dlabel = "Back Down", Xlabel = "Hawkishness",
     cex.lab = .5)+
  ggtitle("Audience Costs (Germany)")

pger_il_hte<-plot(ger_il_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support Torture", 
     Dlabel = "International Law", Xlabel = "Legal Obligation",
     cex.lab = .5)+
  ggtitle("International Law (Germany)")

#India
pind_dp_hte<-plot(ind_dp_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support for Attack", 
     Dlabel = "Democracy", Xlabel = "Support for Democratic Norms",
     cex.lab = .5)+
  ggtitle("Democratic Peace (India)")

pind_ac_hte<-plot(ind_ac_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Leader Approval", 
     Dlabel = "Back Down", Xlabel = "Hawkishness",
     cex.lab = .5)+
  ggtitle("Audience Costs (India)")

pind_il_hte<-plot(ind_il_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support Torture", 
     Dlabel = "International Law", Xlabel = "Legal Obligation",
     cex.lab = .5)+
  ggtitle("International Law (India)")


#Israel
pisr_dp_hte<-plot(isr_dp_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support for Attack", 
     Dlabel = "Democracy", Xlabel = "Support for Democratic Norms",
     cex.lab = .5)+
  ggtitle("Democratic Peace (Israel)")

pisr_il_hte<-plot(isr_il_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support Torture", 
     Dlabel = "International Law", Xlabel = "Legal Obligation",
     cex.lab = .5)+
  ggtitle("International Law (Israel)")

pisr_ac_hte<-plot(isr_ac_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Leader Approval", 
     Dlabel = "Back Down", Xlabel = "Hawkishness",
     cex.lab = .5)+
  ggtitle("Audience Costs (Israel)")

#Japan
pjap_dp_hte<-plot(jap_dp_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support for Attack", 
     Dlabel = "Democracy", Xlabel = "Support for Democratic Norms",
     cex.lab = .5)+
  ggtitle("Democratic Peace (Japan)")

pjap_ac_hte<-plot(jap_ac_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Leader Approval", 
     Dlabel = "Back Down", Xlabel = "Hawkishness",
     cex.lab = .5)+
  ggtitle("Audience Costs (Japan)")

pjap_il_hte<-plot(jap_il_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support Torture", 
     Dlabel = "International Law", Xlabel = "Legal Obligation",
     cex.lab = .5)+
  ggtitle("International Law (Japan)")

#Nigeria
pnig_dp_hte<-plot(nig_dp_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support for Attack", 
     Dlabel = "Democracy", Xlabel = "Support for Democratic Norms",
     cex.lab = .5) + 
  ggtitle("Democratic Peace (Nigeria)")

pnig_ac_hte<-plot(nig_ac_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Leader Approval", 
     Dlabel = "Back Down", Xlabel = "Hawkishness",
     cex.lab = .5)+
  ggtitle("Audience Costs (Nigeria)")

pnig_il_hte<-plot(nig_il_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support Torture", 
     Dlabel = "International Law", Xlabel = "Legal Obligation",
     cex.lab = .5)+
  ggtitle("International Law (Nigeria)")

#USA
pus_dp_hte<-plot(us_dp_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support for Attack", 
     Dlabel = "Democracy", Xlabel = "Support for Democratic Norms",
     cex.lab = .5)+
  ggtitle("Democratic Peace (USA)")

pus_ac_hte<-plot(us_ac_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Leader Approval", 
     Dlabel = "Back Down", Xlabel = "Hawkishness",
     cex.lab = .5)+
ggtitle("Audience Costs (USA)")


pus_il_hte<-plot(us_il_hte,
     theme.bw = T, show.grid = F,
     Ylabel = "Support Torture", 
     Dlabel = "International Law", Xlabel = "Legal Obligation",
     cex.lab = .5)+
ggtitle("International Law (USA)")


(pbr_dp_hte + pbr_ac_hte + pbr_il_hte) /
  (pger_dp_hte + pger_ac_hte + pger_il_hte) /
  (pind_dp_hte + pind_ac_hte + pind_il_hte) /
  (pisr_dp_hte + pisr_ac_hte + pisr_il_hte) /
  (pjap_dp_hte + pjap_ac_hte + pjap_il_hte) / 
  (pnig_dp_hte + pnig_ac_hte + pnig_il_hte) / 
  (pus_dp_hte + pus_ac_hte + pus_il_hte)
ggsave("figures/appendix/hte_bycount.pdf",width = 13, height = 22)


#now report main results for India with Treat*English
eng_ind_dp <- 
  lm(outcome1_dem_new~DP_treat*english,
     data = india_new) 
summary(eng_ind_dp)

eng_ind_ac <- 
  lm(outcome_aud_new~AC_treat*english,
     data = india_new) 
summary(eng_ind_ac)

eng_ind_il <- 
  lm(outcome_il_new~IL_treat*english,
     data = india_new) 
summary(eng_ind_il)

eng_ind_ipe <- 
  lm(outcome_ipe_new~IPE_treat*english,
     data = india_new) 
summary(eng_ind_ipe)

ind_english<-texreg(list(eng_ind_dp,eng_ind_ac,eng_ind_il,eng_ind_ipe),
  label = "tab:ind_english",
  caption.above = T,
  no.margin = T,
  include.ci = FALSE,
  custom.model.names = c("Dem Peace", "Audience Costs", "Int Law","Reciprocity"),
  custom.coef.map = list("IL_treat" = "Intl Law", "english" = "English", 
                         "IL_treat:english" = "IL*English",
                         "AC_treat" = "Back Down",
                         "AC_treat:english" = "BD*English",
                         "DP_treat" = "Democracy",  
                         "DP_treat:english" = "Dem*English",
                         "IPE_treat" = "Harder barrier",  
                         "IPE_treat:english" = "Hard*English"),
  stars = 0.05,
  digits = 3,
  include.rsquared = F, include.rmse = F,
  caption = "Treatment Effect*English in India Sample",
  font.size = "tiny",
  table = T)

ind_english <- gsub("\\begin{tabular}"
               ,"\\resizebox{0.6\\textwidth}{!}{\n\\begin{tabular}"
               ,ind_english
               ,fixed=T)

ind_english <- gsub("\\end{tabular}"
               ,"\\end{tabular}}"
               ,ind_english
               ,fixed=T)

write_file(ind_english, 'tables/ind_english.tex')

# Probing Generalizability ------------

#why do we see strong generalizability? 
#(1) Samples are similar -- let's test for country differences across moderators and covariates

hawk<- lm(hawk~relevel(country_fe, ref = "USA"), all_countries_clean)
legal<- lm(legal~relevel(country_fe, ref = "USA"), all_countries_clean)
norms<- lm(dem_norm~relevel(country_fe, ref = "USA"), all_countries_clean)
age<- lm(age~relevel(country_fe, ref = "USA"), all_countries_clean)
ideology<- lm(as.numeric(ideology)~relevel(country_fe, ref = "USA"), all_countries_clean)
university<- lm(university~relevel(country_fe, ref = "USA"), all_countries_clean)

tab<-stargazer(hawk,legal,norms,age,ideology,university,
                out= "tables/dif_samples.tex",  style="qje",
                title="Estimating Differences Between Country Samples. Each model regresses relevant outcomes over country indicators compared to the US (which serves as a reference category).",
                keep = c("Brazil", "Germany",
                         "India", "Israel", "Japan", "Nigeria"),
                dep.var.labels = c("Hawkishness",
                                   "Legal Oblig",
                                   "Demo Norms",
                                   "Age",
                                   "Ideology",
                                   "University Educated"),
                covariate.labels = c("Brazil", "Germany",
                                     "India", "Israel", "Japan", "Nigeria"),
                omit.stat = c("rsq", "f", "ser", "adj.rsq"),
                star.cutoffs = NULL,
                label = "tab:dif_samples",
                # star.char = c("", "", ""),
                ci=FALSE,
                omit.table.layout = NULL,
                notes.append = FALSE, 
                notes.align = "l",
                notes = c(""))
write(tab, 'tables/dif_samples.tex')


#(2) Limited treatment effect heterogeneity
#include test for systematic heterogeneity, add all pre-treatment variables in the interaction formula

#DP 
HTE_DPBrazil<-all_countries_clean%>%
  filter(country=="Brazil")%>%
  drop_na(outcome1_dem_scaled)%>%
  estimate_systematic(
    outcome1_dem_scaled~DP_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology + university + fear_1+ fear_2,na.rm=T,
    data = .)

summary(HTE_DPBrazil)

HTE_DPGermany<-all_countries_clean%>%
  filter(country=="Germany")%>%
  drop_na(outcome1_dem_scaled)%>%
  estimate_systematic(
    outcome1_dem_scaled~DP_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2,na.rm=T,
    data = .)

summary(HTE_DPGermany)


HTE_DPIndia<-all_countries_clean%>%
  filter(country=="India")%>%
  drop_na(outcome1_dem_scaled)%>%
  estimate_systematic(
    outcome1_dem_scaled~DP_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2,na.rm=T,
    data = .)

summary(HTE_DPIndia)

HTE_DPIsrael<-all_countries_clean%>%
  filter(country=="Israel")%>%
  drop_na(outcome1_dem_scaled)%>%
  estimate_systematic(
    outcome1_dem_scaled~DP_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2,na.rm=T,
    data = .)

summary(HTE_DPIsrael)

HTE_DPJapan<-all_countries_clean%>%
  filter(country=="Japan")%>%
  drop_na(outcome1_dem_scaled)%>%
  estimate_systematic(
    outcome1_dem_scaled~DP_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2,na.rm=T,
    data = .)

summary(HTE_DPJapan)

HTE_DPNigeria<-all_countries_clean%>%
  filter(country=="Nigeria")%>%
  drop_na(outcome1_dem_scaled)%>%
  estimate_systematic(
    outcome1_dem_scaled~DP_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2,na.rm=T,
    data = .)

summary(HTE_DPNigeria)

HTE_DPUSA<-all_countries_clean%>%
  filter(country=="USA")%>%
  drop_na(outcome1_dem_scaled)%>%
  estimate_systematic(
    outcome1_dem_scaled~DP_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2,na.rm=T,
    data = .)

summary(HTE_DPUSA)

#in 3 tests

p.value <- c(HTE_DPUSA$p.value, HTE_DPNigeria$p.value, HTE_DPJapan$p.value, HTE_DPIndia$p.value, 
             HTE_DPIsrael$p.value,HTE_DPGermany$p.value,HTE_DPBrazil$p.value)
Country <- c("USA", "Nigeria", "Japan", "India", "Israel","Germany","Brazil")

df_dp <- data.frame(p.value, Country)

df_dp <-
  df_dp %>% 
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))

df_dp$N_p=ifelse(df_dp$p.value<0.05,1,0)#3
df_dp$N_adjusted=ifelse(df_dp$adjusted<0.05,1,0) #0 after adjusted


#AC 
HTE_ACBrazil<-all_countries_clean%>%
  filter(country=="Brazil")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~IL_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology + university+ fear_1+ fear_2, na.rm=T, 
    data = .)

summary(HTE_ACBrazil)

HTE_ACGermany<-all_countries_clean%>%
  filter(country=="Germany")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~AC_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T, 
    data = .)

summary(HTE_ACGermany)


HTE_ACIndia<-all_countries_clean%>%
  filter(country=="India")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~AC_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology + university+ fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_ACIndia)

HTE_ACIsrael<-all_countries_clean%>%
  filter(country=="Israel")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~AC_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology + university+ fear_1+ fear_2, na.rm=T ,
    data = .)

summary(HTE_ACIsrael)

HTE_ACJapan<-all_countries_clean%>%
  filter(country=="Japan")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~AC_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_ACJapan)

HTE_ACNigeria<-all_countries_clean%>%
  filter(country=="Nigeria")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~AC_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_ACNigeria)

HTE_ACUSA<-all_countries_clean%>%
  filter(country=="USA")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~AC_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_ACUSA)

#3 sig results (Germany, Israel, USA)

p.value <- c(HTE_ACUSA$p.value, HTE_ACNigeria$p.value, HTE_ACJapan$p.value, HTE_ACIndia$p.value, 
             HTE_ACIsrael$p.value,HTE_ACGermany$p.value,HTE_ACBrazil$p.value)

df_ac <- data.frame(p.value, Country)

df_ac <-
  df_ac %>% 
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))

df_ac$N_p=ifelse(df_ac$p.value<0.05,1,0)#3
df_ac$N_adjusted=ifelse(df_ac$adjusted<0.05,1,0) #2 adjusted (IL, GER)


#IL 
HTE_ILBrazil<-all_countries_clean%>%
  filter(country=="Brazil")%>%
  drop_na(outcome_il_scaled)%>%
  estimate_systematic(
    outcome_il_scaled~IL_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_ILBrazil)

HTE_ILGermany<-all_countries_clean%>%
  filter(country=="Germany")%>%
  drop_na(outcome_il_scaled)%>%
  estimate_systematic(
    outcome_il_scaled~IL_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T, 
    data = .)

summary(HTE_ILGermany)


HTE_ILIndia<-all_countries_clean%>%
  filter(country=="India")%>%
  drop_na(outcome_il_scaled)%>%
  estimate_systematic(
    outcome_il_scaled~IL_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_ILIndia)

HTE_ILIsrael<-all_countries_clean%>%
  filter(country=="Israel")%>%
  drop_na(outcome_il_scaled)%>%
  estimate_systematic(
    outcome_il_scaled~IL_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T ,
    data = .)

summary(HTE_ILIsrael)

HTE_ILJapan<-all_countries_clean%>%
  filter(country=="Japan")%>%
  drop_na(outcome_il_scaled)%>%
  estimate_systematic(
    outcome_il_scaled~IL_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology + university+ fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_ILJapan)

HTE_ILNigeria<-all_countries_clean%>%
  filter(country=="Nigeria")%>%
  drop_na(outcome_il_scaled)%>%
  estimate_systematic(
    outcome_il_scaled~IL_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_ILNigeria)

HTE_ILUSA<-all_countries_clean%>%
  filter(country=="USA")%>%
  drop_na(outcome_il_scaled)%>%
  estimate_systematic(
    outcome_il_scaled~IL_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_ILUSA)

#3 sig (USA, Israel, Germany)

p.value <- c(HTE_ILUSA$p.value, HTE_ILNigeria$p.value, HTE_ILJapan$p.value, HTE_ILIndia$p.value, 
             HTE_ILIsrael$p.value,HTE_ILGermany$p.value,HTE_ILBrazil$p.value)

df_il <- data.frame(p.value, Country)

df_il <-
  df_il %>% 
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))

df_il$N_p=ifelse(df_il$p.value<0.05,1,0)#3
df_il$N_adjusted=ifelse(df_il$adjusted<0.05,1,0) #3 

#IPE
HTE_IPEBrazil<-all_countries_clean%>%
  filter(country=="Brazil")%>%
  drop_na(outcome_ipe_scaled)%>%
  estimate_systematic(
    outcome_ipe_scaled~IPE_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_IPEBrazil)

HTE_IPEGermany<-all_countries_clean%>%
  filter(country=="Germany")%>%
  drop_na(outcome_ipe_scaled)%>%
  estimate_systematic(
    outcome_ipe_scaled~IPE_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T, 
    data = .)

summary(HTE_IPEGermany)


HTE_IPEIndia<-all_countries_clean%>%
  filter(country=="India")%>%
  drop_na(outcome_ipe_scaled)%>%
  estimate_systematic(
    outcome_ipe_scaled~IPE_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_IPEIndia)

HTE_IPEIsrael<-all_countries_clean%>%
  filter(country=="Israel")%>%
  drop_na(outcome_ipe_scaled)%>%
  estimate_systematic(
    outcome_ipe_scaled~IPE_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology + university+ fear_1+ fear_2, na.rm=T ,
    data = .)

summary(HTE_IPEIsrael)

HTE_IPEJapan<-all_countries_clean%>%
  filter(country=="Japan")%>%
  drop_na(outcome_ipe_scaled)%>%
  estimate_systematic(
    outcome_ipe_scaled~IPE_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_IPEJapan)

HTE_IPENigeria<-all_countries_clean%>%
  filter(country=="Nigeria")%>%
  drop_na(outcome_ipe_scaled)%>%
  estimate_systematic(
    outcome_ipe_scaled~IPE_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_IPENigeria)

HTE_IPEUSA<-all_countries_clean%>%
  filter(country=="USA")%>%
  drop_na(outcome_ipe_scaled)%>%
  estimate_systematic(
    outcome_ipe_scaled~IPE_treat,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_IPEUSA)

#4 sig

p.value <- c(HTE_IPEUSA$p.value, HTE_IPENigeria$p.value, HTE_IPEJapan$p.value, HTE_IPEIndia$p.value, 
             HTE_IPEIsrael$p.value,HTE_IPEGermany$p.value,HTE_IPEBrazil$p.value)

df_ipe <- data.frame(p.value, Country)

df_ipe <-
  df_ipe %>% 
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))

df_ipe$N_p=ifelse(df_ipe$p.value<0.05,1,0)#4
df_ipe$N_adjusted=ifelse(df_ipe$adjusted<0.05,1,0)#4


#Belligerence cost
HTE_belBrazil<-all_countries_clean%>%
  filter(country=="Brazil")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~belligerence,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T, 
    data = .)

summary(HTE_belBrazil)

HTE_belGermany<-all_countries_clean%>%
  filter(country=="Germany")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~belligerence,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_belGermany)

HTE_belIndia<-all_countries_clean%>%
  filter(country=="India")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~belligerence,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_belIndia)

HTE_belIsrael<-all_countries_clean%>%
  filter(country=="Israel")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~belligerence,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T,
    data = .)

summary(HTE_belIsrael)         

HTE_belJapan<-all_countries_clean%>%
  filter(country=="Japan")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~belligerence,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T, 
    data = .)

summary(HTE_belJapan) 

HTE_belNigeria<-all_countries_clean%>%
  filter(country=="Nigeria")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~belligerence,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T, 
    data = .)

summary(HTE_belNigeria) 

HTE_belUSA<-all_countries_clean%>%
  filter(country=="USA")%>%
  drop_na(outcome_aud_scaled)%>%
  estimate_systematic(
    outcome_aud_scaled~belligerence,
    interaction.formula = ~ gender_fe + age + dem_norm + hawk + legal +
      voting_fe + ideology+ university + fear_1+ fear_2, na.rm=T, 
    data = .)

summary(HTE_belUSA) 

#7 sig tests (all of them..)

p.value <- c(HTE_belUSA$p.value, HTE_belNigeria$p.value, HTE_belJapan$p.value, HTE_belIndia$p.value, 
             HTE_belIsrael$p.value,HTE_belGermany$p.value,HTE_belBrazil$p.value)

df_bel <- data.frame(p.value, Country)

df_bel <-
  df_bel %>% 
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))

df_bel$N_p=ifelse(df_bel$p.value<0.05,1,0)#7
df_bel$N_adjusted=ifelse(df_bel$adjusted<0.05,1,0)#7

##Descriptive Statistics------

all_countries_clean$info_leak_ipe<-ifelse(all_countries_clean$info_leak_ipe=="",NA,
                                          all_countries_clean$info_leak_ipe)

all_attentive<-all_countries_clean%>%
  filter(info_leak_ipe!="")%>%
  mutate(.,
         vote = case_when(voting_notWVS=="1"~1,
                           voting_notWVS=="2"~0),
         gender = case_when(S=="1"~1,
                            S=="2"~0),
         education_new = as.numeric(education))


disc_table_all <- dplyr::select(all_attentive, 
                            c(outcome1_dem_new,outcome2_dem_new,
                              outcome_aud_new,outcome_il_new,
                              outcome_ipe_new,DP_man,AC_man,IL_man,IPE_man,
                              dem_norm,hawk,legal,gender, education_new,
                              vote,age))

stargazer(as.data.frame(disc_table_all),
          covariate.labels = c("DP outcome","DP outcome 2",
                               "AC outcome","IL outcome",
                               "FDI outcome","Manipulation DP",
                               "Manipulation AC","Manipulation IL","Manipulation FDI",
                               "Democratic norms","Hawkishness","Intl legal obligation",
                               "Gender","Education","Eligable to vote","Age"),
          
          title = "Descriptive Statistics - All Countries",
          label = "tab:dsc_all",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          out="tables/desc_all.tex")


#create table for each country
#Brazil
brazil_attentive<-all_attentive%>%
  filter(country=="Brazil")

disc_table_brazil <- dplyr::select(brazil_attentive, 
                                  c(outcome1_dem_new,outcome2_dem_new,
                                    outcome_aud_new,outcome_il_new,
                                    outcome_ipe_new,DP_man,AC_man,IL_man,IPE_man,
                                    dem_norm,hawk,legal,gender,education_new,
                                    vote,age))

stargazer(as.data.frame(disc_table_brazil),
          covariate.labels = c("DP outcome","DP outcome 2",
                               "AC outcome","IL outcome",
                               "FDI outcome","Manipulation DP",
                               "Manipulation AC","Manipulation IL","Manipulation FDI",
                               "Democratic norms","Hawkishness","Intl legal obligation",
                               "Gender","Education","Eligable to vote","Age"),
          
          title = "Descriptive Statistics - Brazil",
          label = "tab:dsc_br",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          out="tables/desc_br.tex")

#Germany
germany_attentive<-all_attentive%>%
  filter(country=="Germany")

disc_table_germany <- dplyr::select(germany_attentive, 
                                    c(outcome1_dem_new,outcome2_dem_new,
                                      outcome_aud_new,outcome_il_new,
                                      outcome_ipe_new,DP_man,AC_man,IL_man,IPE_man,
                                      dem_norm,hawk,legal,gender,education_new,
                                      vote,age))

stargazer(as.data.frame(disc_table_germany),
          covariate.labels = c("DP outcome","DP outcome 2",
                               "AC outcome","IL outcome",
                               "FDI outcome","Manipulation DP",
                               "Manipulation AC","Manipulation IL","Manipulation FDI",
                               "Democratic norms","Hawkishness","Intl legal obligation",
                               "Gender","Education","Eligable to vote","Age"),
          
          title = "Descriptive Statistics - Germany",
          label = "tab:dsc_gr",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          out="tables/desc_gr.tex")
  
#India
india_attentive<-all_attentive%>%
  filter(country=="India")

disc_table_india <- dplyr::select(india_attentive, 
                                    c(outcome1_dem_new,outcome2_dem_new,
                                      outcome_aud_new,outcome_il_new,
                                      outcome_ipe_new,DP_man,AC_man,IL_man,IPE_man,
                                      dem_norm,hawk,legal,gender,education_new,
                                      vote,age))

stargazer(as.data.frame(disc_table_india),
          covariate.labels = c("DP outcome","DP outcome 2",
                               "AC outcome","IL outcome",
                               "FDI outcome","Manipulation DP",
                               "Manipulation AC","Manipulation IL","Manipulation FDI",
                               "Democratic norms","Hawkishness","Intl legal obligation",
                               "Gender","Education","Eligable to vote","Age"),
          
          title = "Descriptive Statistics - India",
          label = "tab:dsc_in",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          out="tables/desc_in.tex")  

#Israel
israel_attentive<-all_attentive%>%
  filter(country=="Israel")

disc_table_israel <- dplyr::select(israel_attentive, 
                                  c(outcome1_dem_new,outcome2_dem_new,
                                    outcome_aud_new,outcome_il_new,
                                    outcome_ipe_new,DP_man,AC_man,IL_man,IPE_man,
                                    dem_norm,hawk,legal,gender,education_new,
                                    vote,age))

stargazer(as.data.frame(disc_table_israel),
          covariate.labels = c("DP outcome","DP outcome 2",
                               "AC outcome","IL outcome",
                               "FDI outcome","Manipulation DP",
                               "Manipulation AC","Manipulation IL","Manipulation FDI",
                               "Democratic norms","Hawkishness","Intl legal obligation",
                               "Gender","Education","Eligable to vote","Age"),
          
          title = "Descriptive Statistics - Israel",
          label = "tab:dsc_il",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          out="tables/desc_il.tex")  


#Japan
japan_attentive<-all_attentive%>%
  filter(country=="Japan")

disc_table_japan <- dplyr::select(japan_attentive, 
                                   c(outcome1_dem_new,outcome2_dem_new,
                                     outcome_aud_new,outcome_il_new,
                                     outcome_ipe_new,DP_man,AC_man,IL_man,IPE_man,
                                     dem_norm,hawk,legal,gender,education_new,
                                     vote,age))

stargazer(as.data.frame(disc_table_japan),
          covariate.labels = c("DP outcome","DP outcome 2",
                               "AC outcome","IL outcome",
                               "FDI outcome","Manipulation DP",
                               "Manipulation AC","Manipulation IL","Manipulation FDI",
                               "Democratic norms","Hawkishness","Intl legal obligation",
                               "Gender","Education","Eligable to vote","Age"),
          
          title = "Descriptive Statistics - Japan",
          label = "tab:dsc_jp",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          out="tables/desc_jp.tex")  


#Nigeria
nigeria_attentive<-all_attentive%>%
  filter(country=="Nigeria")

disc_table_nigeria <- dplyr::select(nigeria_attentive, 
                                  c(outcome1_dem_new,outcome2_dem_new,
                                    outcome_aud_new,outcome_il_new,
                                    outcome_ipe_new,DP_man,AC_man,IL_man,IPE_man,
                                    dem_norm,hawk,legal,gender,education_new,
                                    vote,age))

stargazer(as.data.frame(disc_table_nigeria),
          covariate.labels = c("DP outcome","DP outcome 2",
                               "AC outcome","IL outcome",
                               "FDI outcome","Manipulation DP",
                               "Manipulation AC","Manipulation IL","Manipulation FDI",
                               "Democratic norms","Hawkishness","Intl legal obligation",
                               "Gender","Education","Eligable to vote","Age"),
          
          title = "Descriptive Statistics - Nigeria",
          label = "tab:dsc_ng",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          out="tables/desc_ng.tex") 

#USA
us_attentive<-all_attentive%>%
  filter(country=="USA")

disc_table_us <- dplyr::select(us_attentive, 
                                    c(outcome1_dem_new,outcome2_dem_new,
                                      outcome_aud_new,outcome_il_new,
                                      outcome_ipe_new,DP_man,AC_man,IL_man,IPE_man,
                                      dem_norm,hawk,legal,gender,education_new,
                                      vote,age))

stargazer(as.data.frame(disc_table_us),
          covariate.labels = c("DP outcome","DP outcome 2",
                               "AC outcome","IL outcome",
                               "FDI outcome","Manipulation DP",
                               "Manipulation AC","Manipulation IL","Manipulation FDI",
                               "Democratic norms","Hawkishness","Intl legal obligation",
                               "Gender","Education","Eligable to vote","Age"),
          
          title = "Descriptive Statistics - USA",
          label = "tab:dsc_us",
          style = "qje",
          notes.append = T,
          notes.align = "l",
          out="tables/desc_us.tex") 


##DP extention ------
#We examine the effect of democracy on a second outcome -- support joining an intl mission
us_dp2 <- lm_robust(outcome2_dem_new ~ DP_treat, subset = country == "USA", data = all_countries_clean)
br_dp2 <- lm_robust(outcome2_dem_new ~ DP_treat, subset = country == "Brazil", data = all_countries_clean)
gr_dp2 <- lm_robust(outcome2_dem_new ~ DP_treat, subset = country == "Germany", data = all_countries_clean)
in_dp2 <- lm_robust(outcome2_dem_new ~ DP_treat, subset = country == "India", data = all_countries_clean)
il_dp2 <- lm_robust(outcome2_dem_new ~ DP_treat, subset = country == "Israel", data = all_countries_clean)
jp_dp2 <- lm_robust(outcome2_dem_new ~ DP_treat, subset = country == "Japan", data = all_countries_clean)
ng_dp2 <- lm_robust(outcome2_dem_new ~ DP_treat, subset = country == "Nigeria", data = all_countries_clean)

#print table
texreg(list(
  br_dp2,gr_dp2,in_dp2,il_dp2,jp_dp2,ng_dp2,us_dp2),
  label = "tab:treat_outcome2",
  caption.above = T,
  no.margin = T,
  include.ci = FALSE,
  custom.header = list("Joint International Mission" = 1:7),
  custom.model.names	= c("BRZ", "GRM", "IND", "ISL", "JPN", "NGR", "USA"),
  custom.coef.map = list("DP_treat" = "Democracy"),
  stars = 0.05,
  digits = 2,
  include.rsquared = F, include.rmse = F,
  caption = "Alternative Outcome (DP): Treatment Effect on Joining Intl Mission",
  font.size = "scriptsize",
  table = T,
  file = "tables/outcome2_DP.tex")


#AC extension ------
#We break down belligerence cost and inconsistency cost

#belligerence cost
BEL_bycountry <-
  all_countries_clean %>% 
  group_by(country) %>%
  do(tidy(lm_robust(outcome_aud_new ~ belligerence, data = .)))%>%
  filter(term!="(Intercept)")%>%
  mutate(.,
         meta = "Country Specific Replication",
         study = "Belligerence Costs")%>%ungroup()%>%
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))


##Estimate Meta-Analytic Random Effects Model:
BEL_meta <- 
  rma(yi=estimate, sei=std.error, data= BEL_bycountry,
      slab = country) %>% tidy(.,  conf.int=T) %>% 
  mutate(.,
         country = "All Countries", 
         meta = "Meta",
         study = "Belligerence Costs",
         df = sum(BEL_bycountry$df))

#BEL_meta$I2 #98.5

#inconsistency cost
INC_bycountry <-
  all_countries_clean %>% 
  group_by(country) %>%
  do(tidy(lm_robust(outcome_aud_new ~ inconsistency, data = .)))%>%
  filter(term!="(Intercept)")%>%
  mutate(.,
         meta = "Country Specific Replication",
         study = "Inconsistency Costs")%>%ungroup()%>%
  mutate(.,
         adjusted = p.adjust(p.value, method = "BH"))


##Estimate Meta-Analytic Random Effects Model:
INC_meta <- 
  rma(yi=estimate, sei=std.error, data= INC_bycountry,
      slab = country) %>% tidy(.,  conf.int=T) %>% 
  mutate(.,
         country = "All Countries", 
         meta = "Meta",
         study = "Inconsistency Costs",
         df = sum(INC_bycountry$df))

#Estimate ATE from BK original study:
#`Strategy` is treatment where (3 not_engage), (2 engage), (1 stay_out)
audience_orig <- 
  audience_orig %>%
  mutate(.,
         bel = case_when(
           Strategy == 2 ~ 1, 
           Strategy == 1 ~ 0),
         inconsist = case_when(
           Strategy == 3 ~ 1, 
           Strategy == 2 ~ 0))

# Estimate original models
BEL_orig <-
  audience_orig %>% 
  lm_robust(totApprov1~
              bel,
            data = .) %>% 
  tidy(.) %>% 
  filter(term == "bel") %>% 
  mutate(.,
         meta = "Original",
         country = "USA",
         study = "Belligerence Costs")

INC_orig <-
  audience_orig %>% 
  lm_robust(totApprov1~
              inconsist,
            data = .) %>% 
  tidy(.) %>% 
  filter(term == "inconsist") %>% 
  mutate(.,
         meta = "Original",
         country = "USA",
         study = "Inconsistency Costs")

#combine models
AC_extensions <-
  bind_rows(BEL_bycountry,BEL_meta,BEL_orig,
            INC_bycountry,INC_meta,INC_orig,)%>%
  mutate(.,
         p_report = ifelse(meta=="Country Specific Replication",adjusted,p.value))

AC_extensions$country <- 
  factor(AC_extensions$country,
         levels = c("USA","Nigeria","Japan", "India",
                    "Israel","Germany", "Brazil", "All Countries"))

AC_extensions$meta <- 
  factor(AC_extensions$meta,
         levels = c("Meta","Country Specific Replication","Original"))

ggplot(AC_extensions, aes(x = as.numeric(estimate), y =as.factor(country), 
                               color =meta, shape = meta)) +
  geom_vline(xintercept = 0, color = "dodgerblue4", linetype = 2, size = 0.2) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high)) +
  geom_text(aes(label= paste("p =", format(round(p_report, digits = 2), nsmall = 2))), hjust=0.3, vjust=-1.4, size = 3.3) +
  #geom_text(aes(label= paste("p. adj = ",format(round(adjusted, digits = 2), nsmall =2))), hjust=-.33, vjust=-0.5, size = 3.3, 
  #         family = "Times")+
  labs(x = "Main Effect",
       y = "") +
  facet_grid(cols = vars(study), rows = vars(meta), scales="free_y", space = "free")+
  scale_color_manual(values = c("grey40", "firebrick4", "dodgerblue4"))+
  xlim(-1.1, 1.1)+
  theme_bw()+
  theme(strip.background =element_rect(fill="dodgerblue4"),
        strip.text = element_text(colour = 'white'),
        panel.grid.major = element_blank(),
        legend.position = "none",
        strip.text.x = element_text(size = 12),
        strip.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        axis.text.x = element_text(size = 11),
        axis.text.y = element_text(size = 11),
        panel.grid.minor = element_blank(),
        panel.background = element_blank())
ggsave("figures/appendix/AC_extension.pdf", width = 6, height = 6)


#create table for appendix 
AC_extensions$country_n<-ifelse(AC_extensions$meta=="Original",
                                     "Original (USA)",as.character(AC_extensions$country))

SE_beli<-subset(AC_extensions$std.error,AC_extensions$study=="Belligerence Costs")
SE_incon<-subset(AC_extensions$std.error,AC_extensions$study=="Inconsistency Costs")

Estimate_beli<-subset(AC_extensions$estimate,AC_extensions$study=="Belligerence Costs")
Estimate_incon<-subset(AC_extensions$estimate,AC_extensions$study=="Inconsistency Costs")

P_beli<-subset(AC_extensions$p_report,AC_extensions$study=="Belligerence Costs")
P_incon<-subset(AC_extensions$p_report,AC_extensions$study=="Inconsistency Costs")

N_beli<-subset(AC_extensions$df,AC_extensions$study=="Belligerence Costs")
N_incon<-subset(AC_extensions$df,AC_extensions$study=="Inconsistency Costs")


AC_extensions_wide<-data.frame(Samples = unique(AC_extensions$country_n),
                               Estimate_beli = Estimate_beli,SE_beli = SE_beli,P_beli = P_beli,N_beli = N_beli,
                               Estimate_incon=Estimate_incon,SE_incon=SE_incon,P_incon=P_incon,N_incon=N_incon)

AC_extensions_wide<-rename(AC_extensions_wide,
                         "Estimate (Belligerence)" = "Estimate_beli",
                         "Estimate (Inconsistency)" = "Estimate_incon",
                         
                         "SE (Belligerence)" = "SE_beli",
                         "SE (Inconsistency)" = "SE_incon",
                         
                         "P value (Belligerence)" = "P_beli",
                         "P value (Inconsistency)" = "P_incon",
                         
                         "N (DP)" = "N_beli",
                         "N (Inconsistency)" = "N_incon")

write(kable(AC_extensions_wide, format = "latex", booktabs = T, caption = 
              "Audience Costs Extension (Figure A8) in table form. \\label{tab:ACextable}", digits=2) %>% 
        add_header_above(c(" ", "Beligerence Costs" = 4, "Inconsistency Costs" = 4)) %>% 
        kable_styling(latex_options="scale_down"),'tables/ACextable.tex')


#examine the moderation of hawkishness by belligerence and inconsistency costs

#HTE inconsistency and hawkishness
all_countries_clean<-as.data.frame(all_countries_clean)

mod_inc<- interflex(Y = "outcome_aud_new", D = "inconsistency", X = "hawk", 
                      Z = c("country_fe", "gender_fe", "age","education_fe", "voting_fe"), data = all_countries_clean,
                      estimator = "binning", vcov.type = "robust",
                      na.rm = T,
                      main = "Marginal Effects")


plot(mod_inc,
     theme.bw = T, show.grid = F,
     Ylabel = "Leader Approval", 
     Dlabel = "Inconsistency", Xlabel = "Hawkishness",
     cex.lab = 1)

ggsave("figures/appendix/mod_inc.pdf",width = 6, height = 6)


#HTE belligerence and hawkishness
mod_bel<- interflex(Y = "outcome_aud_new", D = "belligerence", X = "hawk", 
                    Z = c("country_fe", "gender_fe", "age","education_fe", "voting_fe"), data = all_countries_clean,
                    estimator = "binning", vcov.type = "robust",
                    na.rm = T,
                    main = "Marginal Effects")

plot(mod_bel,
     theme.bw = T, show.grid = F,
     Ylabel = "Leader Approval", 
     Dlabel = "Belligerence", Xlabel = "Hawkishness",
     cex.lab = 1)

ggsave("figures/appendix/mod_bel.pdf",width = 6, height = 6)

#create appendix-table for mod_bel.pdf

bel_HTE<-lm(outcome_aud_new~belligerence*hawk+country_fe+gender_fe+age+education_fe+voting_fe,all_countries_clean)
summary(bel_HTE)

inc_HTE<-lm(outcome_aud_new~inconsistency*hawk+country_fe+gender_fe+age+education_fe+voting_fe,all_countries_clean)
summary(inc_HTE)

texreg(list(
  bel_HTE,inc_HTE),
  label = "tab:mod_bel",
  caption.above = T,
  no.margin = T,
  include.ci = FALSE,
  custom.model.names	= c("Belligerence Costs","Inconsistency Costs" ),
  custom.coef.map = list("belligerence:hawk" = "Belligerence*Hawk",
                         "belligerence" = "Belligerence",
                         "hawk"="Hawk", "inconsistency:hawk"="Inconsistency*Hawk",
                         "inconsistency"="Inconsistency"),
  stars = 0.05,
  digits = 2,
  include.rsquared = F, include.rmse = F,
  caption = "Moderating effect of Hawkishness in AC extension (Figure A9) in table form.",
  font.size = "scriptsize",
  table = T,
  file = "tables/mod_belTab.tex")

#repeat belligerence cost by country

#USA
us_bel<-interflex(Y = "outcome_aud_new", D = "belligerence", X = "hawk", 
                            Z = c("gender_fe", "age","education_fe", "voting_fe"), data = us_new,
                            estimator = "binning", vcov.type = "robust",
                            na.rm = T,
                            main = "Marginal Effects")

#Nigeria
nig_bel<- interflex(Y = "outcome_aud_new", D = "belligerence", X = "hawk", 
                    Z = c("gender_fe", "age","education_fe", "voting_fe"), data = nig_new,
                    estimator = "binning", vcov.type = "robust",
                    na.rm = T,
                    main = "Marginal Effects")

#Japan
jap_bel<- interflex(Y = "outcome_aud_new", D = "belligerence", X = "hawk", 
                    Z = c("gender_fe", "age","education_fe", "voting_fe"), data = jap_new,
                    estimator = "binning", vcov.type = "robust",
                    na.rm = T,
                    main = "Marginal Effects")

#Israel
isr_bel<- interflex(Y = "outcome_aud_new", D = "belligerence", X = "hawk", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = israel_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

#India
ind_bel<- interflex(Y = "outcome_aud_new", D = "belligerence", X = "hawk", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = india_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

#Germany
ger_bel<- interflex(Y = "outcome_aud_new", D = "belligerence", X = "hawk", 
                       Z = c("gender_fe", "age","education_fe", "voting_fe"), data = germany_new,
                       estimator = "binning", vcov.type = "robust",
                       na.rm = T,
                       main = "Marginal Effects")

#Brazil
br_bel<- interflex(Y = "outcome_aud_new", D = "belligerence", X = "hawk", 
                      Z = c("gender_fe", "age","education_fe", "voting_fe"), data = brazil_new,
                      estimator = "binning", vcov.type = "robust",
                      na.rm = T,
                      main = "Marginal Effects")

library(patchwork)

#plot

pus_bel<-plot(us_bel,
                 theme.bw = T, show.grid = F,
                 Ylabel = "Leader Approval", 
                 Dlabel = "Belligerence", Xlabel = "Hawkishness",
                 cex.lab = .5) +
  ggtitle('USA')

pnig_bel<-plot(nig_bel,
              theme.bw = T, show.grid = F,
              Ylabel = "Leader Approval", 
              Dlabel = "Belligerence", Xlabel = "Hawkishness",
              cex.lab = .5) +
  ggtitle('Nigeria')

pjap_bel<-plot(jap_bel,
               theme.bw = T, show.grid = F,
               Ylabel = "Leader Approval", 
               Dlabel = "Belligerence", Xlabel = "Hawkishness",
               cex.lab = .5) +
  ggtitle('Japan')


pisr_bel<-plot(isr_bel,
               theme.bw = T, show.grid = F,
               Ylabel = "Leader Approval", 
               Dlabel = "Belligerence", Xlabel = "Hawkishness",
               cex.lab = .5) +
  ggtitle('Israel')


pind_bel<-plot(ind_bel,
               theme.bw = T, show.grid = F,
               Ylabel = "Leader Approval", 
               Dlabel = "Belligerence", Xlabel = "Hawkishness",
               cex.lab = .5) +
  ggtitle('India')


pger_bel<-plot(ger_bel,
               theme.bw = T, show.grid = F,
               Ylabel = "Leader Approval", 
               Dlabel = "Belligerence", Xlabel = "Hawkishness",
               cex.lab = .5) +
  ggtitle('Germany')


pbr_bel<-plot(br_bel,
               theme.bw = T, show.grid = F,
               Ylabel = "Leader Approval", 
               Dlabel = "Belligerence", Xlabel = "Hawkishness",
               cex.lab = .5) +
  ggtitle('Brazil')


(pus_bel + pnig_bel + pjap_bel) /
  (pisr_bel + pind_bel) / 
  (pger_bel+ pbr_bel)
ggsave("figures/appendix/bel_HTE_bycount.pdf",width = 13, height = 15)

#create table form for bel_HTE_bycount.pdf 

us_beli<-lm(outcome_aud_new~belligerence*hawk+gender_fe+age+education_fe+voting_fe,us_new)
nig_beli<-lm(outcome_aud_new~belligerence*hawk+gender_fe+age+education_fe+voting_fe,nig_new)
jap_beli<-lm(outcome_aud_new~belligerence*hawk+gender_fe+age+education_fe+voting_fe,jap_new)
israel_beli<-lm(outcome_aud_new~belligerence*hawk+gender_fe+age+education_fe+voting_fe,israel_new)
india_beli<-lm(outcome_aud_new~belligerence*hawk+gender_fe+age+education_fe+voting_fe,india_new)
brazil_beli<-lm(outcome_aud_new~belligerence*hawk+gender_fe+age+education_fe+voting_fe,brazil_new)
germany_beli<-lm(outcome_aud_new~belligerence*hawk+gender_fe+age+education_fe+voting_fe,germany_new)


texreg(list(brazil_beli,germany_beli,india_beli,israel_beli,jap_beli,nig_beli,us_beli),
  label = "tab:bel_HTE_bycount",
  caption.above = T,
  no.margin = T,
  include.ci = FALSE,
  custom.model.names	= c("BRZ", "GRM", "IND", "ISL", "JPN", "NGR", "USA"),
  custom.coef.map = list("belligerence:hawk" = "Belligerence*Hawk",
                         "belligerence" = "Belligerence",
                         "hawk"="Hawk", "inconsistency:hawk"="Inconsistency*Hawk",
                         "inconsistency"="Inconsistency"),
  stars = 0.05,
  digits = 2,
  include.rsquared = F, include.rmse = F,
  caption = "Belligerence costs moderating effects by country (Figure A10) in table form.",
  font.size = "scriptsize",
  table = T,
  file = "tables/bel_HTE_bycountTAB.tex")


#Information leakage ------
#are treated respondents more likely to think of a specific country?

#DP
us_dp_info <- lm_robust(DP_info ~ DP_treat, subset = country == "USA", data = all_countries_clean)
br_dp_info <- lm_robust(DP_info ~ DP_treat, subset = country == "Brazil", data = all_countries_clean)
gr_dp_info <- lm_robust(DP_info ~ DP_treat, subset = country == "Germany", data = all_countries_clean)
in_dp_info <- lm_robust(DP_info ~ DP_treat, subset = country == "India", data = all_countries_clean)
il_dp_info <- lm_robust(DP_info ~ DP_treat, subset = country == "Israel", data = all_countries_clean)
jp_dp_info <- lm_robust(DP_info ~ DP_treat, subset = country == "Japan", data = all_countries_clean)
ng_dp_info <- lm_robust(DP_info ~ DP_treat, subset = country == "Nigeria", data = all_countries_clean)

#AC
us_ac_info <- lm_robust(AC_info ~ AC_treat, subset = country == "USA", data = all_countries_clean)
br_ac_info <- lm_robust(AC_info ~ AC_treat, subset = country == "Brazil", data = all_countries_clean)
gr_ac_info <- lm_robust(AC_info ~ AC_treat, subset = country == "Germany", data = all_countries_clean)
in_ac_info <- lm_robust(AC_info ~ AC_treat, subset = country == "India", data = all_countries_clean)
il_ac_info <- lm_robust(AC_info ~ AC_treat, subset = country == "Israel", data = all_countries_clean)
jp_ac_info <- lm_robust(AC_info ~ AC_treat, subset = country == "Japan", data = all_countries_clean)
ng_ac_info <- lm_robust(AC_info ~ AC_treat, subset = country == "Nigeria", data = all_countries_clean)


#IL
us_il_info <- lm_robust(IL_info ~ IL_treat, subset = country == "USA", data = all_countries_clean)
br_il_info <- lm_robust(IL_info ~ IL_treat, subset = country == "Brazil", data = all_countries_clean)
gr_il_info <- lm_robust(IL_info ~ IL_treat, subset = country == "Germany", data = all_countries_clean)
in_il_info <- lm_robust(IL_info ~ IL_treat, subset = country == "India", data = all_countries_clean)
il_il_info <- lm_robust(IL_info ~ IL_treat, subset = country == "Israel", data = all_countries_clean)
jp_il_info <- lm_robust(IL_info ~ IL_treat, subset = country == "Japan", data = all_countries_clean)
ng_il_info <- lm_robust(IL_info ~ IL_treat, subset = country == "Nigeria", data = all_countries_clean)


#IPE
us_ip_info <- lm_robust(IPE_info ~ IPE_treat, subset = country == "USA", data = all_countries_clean)
br_ip_info <- lm_robust(IPE_info ~ IPE_treat, subset = country == "Brazil", data = all_countries_clean)
gr_ip_info <- lm_robust(IPE_info ~ IPE_treat, subset = country == "Germany", data = all_countries_clean)
in_ip_info <- lm_robust(IPE_info ~ IPE_treat, subset = country == "India", data = all_countries_clean)
il_ip_info <- lm_robust(IPE_info ~ IPE_treat, subset = country == "Israel", data = all_countries_clean)
jp_ip_info <- lm_robust(IPE_info ~ IPE_treat, subset = country == "Japan", data = all_countries_clean)
ng_ip_info <- lm_robust(IPE_info ~ IPE_treat, subset = country == "Nigeria", data = all_countries_clean)

#print table
texreg(list(
  us_dp_info,br_dp_info,gr_dp_info,in_dp_info,il_dp_info,jp_dp_info,ng_dp_info,
  us_ac_info,br_ac_info,gr_ac_info,in_ac_info,il_ac_info,jp_ac_info,ng_ac_info,
  us_il_info,br_il_info,gr_il_info,in_il_info,il_il_info,jp_il_info,ng_il_info,
  us_ip_info,br_ip_info,gr_ip_info,in_ip_info,il_ip_info,jp_ip_info,ng_ip_info),
  label = "tab:treat_info",
  caption.above = T,
  no.margin = T,
  include.ci = FALSE,
    custom.header = list("Info leak DP" = 1:7,
                         "Info leak AC" = 8:14,
                         "Info leak IL" = 15:21,
                         "Info leak FDI" = 22:28),
  custom.model.names	= c("USA", "BRZ", "GRM", "IND", "ISL", "JPN", "NGR",
                         "USA", "BRZ", "GRM", "IND", "ISL", "JPN", "NGR",
                         "USA", "BRZ", "GRM", "IND", "ISL", "JPN", "NGR",
                         "USA", "BRZ", "GRM", "IND", "ISL", "JPN", "NGR"),
  custom.coef.map = list("DP_treat" = "Democracy", "AC_treat" = "Engage",  
                         "IL_treat" = "IL Law", "IPE_treat" = "Harder barrier"),
  fontsize = "tiny",
  stars = 0.05,
  digits = 2,
  include.rsquared = F, include.rmse = F,
  caption = "Information Leakage: ATE on thinking of a specific country",
  font.size = "scriptsize",
  table = T,
  file = "tables/info.tex")


#Plot info leakage by experiment----

#DP
info_DP<-all_countries_clean%>%
  dplyr::select(country,DP_info1)%>%
  mutate(.,
         experiment = "Democratic Peace")%>%
  rename(info = DP_info1)%>% drop_na() 
#AC
info_AC<-all_countries_clean%>%
  dplyr::select(country,AC_info1)%>%
  mutate(.,
         experiment = "Audience Costs")%>%
  rename(info = AC_info1)%>% drop_na()
#IL
info_IL<-all_countries_clean%>%
  dplyr::select(country,IL_info1)%>%
  mutate(.,
         experiment = "International Law")%>%
  rename(info = IL_info1)%>% drop_na()
#IPE
info_IPE<-all_countries_clean%>%
  dplyr::select(country,IPE_info1)%>%
  mutate(.,
         experiment = "Receprocity FDI")%>%
  rename(info = IPE_info1)%>% drop_na()

info_all<-rbind(info_DP,info_AC,
                 info_IL,info_IPE)


ggplot(info_all, aes(x= info,  group=country)) + 
  geom_bar(aes(y = ..prop..), stat="count") +
  geom_text(aes(label = scales::percent(..prop..),
                y= ..prop.. ), stat= "count", vjust = -.5, size = 2) +
  labs(y = "Percent", x="Information Leakage") +
  facet_grid(country~experiment) +
  theme_bw()+
  theme(text = element_text(size = 14, family = "Times"),
        legend.position = "none",
        panel.grid.major = element_blank(),
        panel.background = element_blank()
  )+
  scale_y_continuous(labels = scales::percent, limits = c(0,1))

ggsave("figures/appendix/info_leak.pdf",width = 14, height = 7)


#plot top 5 words per condition per country by experiment

#clean texts

all_countries_clean$info_leak_aud_1_TEXT<-change_words(all_countries_clean$info_leak_aud_1_TEXT)
all_countries_clean$info_leak_dem_1_TEXT<-change_words(all_countries_clean$info_leak_dem_1_TEXT)
all_countries_clean$info_leak_il_1_TEXT<-change_words(all_countries_clean$info_leak_il_1_TEXT)
all_countries_clean$info_leak_ipe_1_TEXT<-change_words(all_countries_clean$info_leak_ipe_1_TEXT)

#DP
clean_text_DP <- all_countries_clean %>%
  filter(info_leak_dem_1_TEXT!="") %>%
  dplyr::select(info_leak_dem_1_TEXT,country,dem_condition)%>%
  unnest_tokens(word, info_leak_dem_1_TEXT)%>%
  anti_join(stop_words, by = "word") %>%
  group_by(country,dem_condition) %>%
  count(word) %>% 
  top_n(5, n) %>%
  arrange(desc(n), .by_group = TRUE)%>%
  mutate(.,
         new_word = factor(word, levels=word[order(n,decreasing=F)]))

myplots <- lapply(split(clean_text_DP,list(clean_text_DP$country,clean_text_DP$dem_condition)), function(x){
  x$word <- factor(x$word, levels=x$word[order(x$n,decreasing=F)])
  x$dem_condition <- ifelse(x$dem_condition=="democracy","Democracy","Non Democracy")
  x$country_lab<-ifelse(x$dem_condition=="Democracy",x$country,"")
  #make the plot
  p <- ggplot(x, aes(x = word, y = n, width=0.75)) +
    geom_bar(stat = "identity") + #to force all levels to be considered
    ylim(0,900)+
    theme_bw()+
    theme(legend.position="none")+
    labs(y="count", x=unique(x$country_lab),
         title = unique(x$dem_condition))+
    coord_flip()
})

myplots <- myplots[order(names(myplots))]
a<-do.call(arrangeGrob,(c(myplots, ncol=2)))
ggsave("figures/appendix/words_DP.pdf",a, width = 8, height = 12)


#AC
clean_text_AC <- all_countries_clean %>%
  filter(info_leak_aud_1_TEXT!="") %>%
  dplyr::select(info_leak_aud_1_TEXT,country,aud_condition)%>%
  unnest_tokens(word, info_leak_aud_1_TEXT)%>% 
  anti_join(stop_words, by = "word") %>%
  group_by(country,aud_condition) %>%
  count(word) %>% 
  top_n(5, n) %>%
  arrange(desc(n), .by_group = TRUE)%>%
  mutate(.,
         new_word = factor(word, levels=word[order(n,decreasing=F)]))


myplotsAC <- lapply(split(clean_text_AC,list(clean_text_AC$country,clean_text_AC$aud_condition)), function(x){
  x$word <- factor(x$word, levels=x$word[order(x$n,decreasing=F)])
  x<-x%>%
    mutate(.,
           aud_condition=case_when(
             aud_condition=="engage"~"Engage",
             aud_condition=="not_engage"~"Not Engage",
             aud_condition=="stay_out"~"Stay Out"))
  x$country_lab<-ifelse(x$aud_condition=="Engage",x$country,"")
 
  #make the plot
  p <- ggplot(x, aes(x = word, y = n, width=0.75)) +
    geom_bar(stat = "identity") + #to force all levels to be considered
    theme_bw()+
    ylim(0,500)+
    theme(legend.position="none")+
    labs(y="count", x=unique(x$country_lab),
         title = unique(x$aud_condition))+
    coord_flip()
})

myplotsAC <- myplotsAC[order(names(myplotsAC))]
a<-do.call(arrangeGrob,(c(myplotsAC, ncol=3)))
ggsave("figures/appendix/words_AC.pdf",a, width = 8, height = 10)


#IL
clean_text_IL <- all_countries_clean %>%
  filter(info_leak_il_1_TEXT!="") %>%
  dplyr::select(info_leak_il_1_TEXT,country,il_condition)%>%
  unnest_tokens(word, info_leak_il_1_TEXT)%>% 
  anti_join(stop_words, by = "word") %>%
  group_by(country,il_condition) %>%
  count(word) %>% 
  top_n(5, n) %>%
  arrange(desc(n), .by_group = TRUE)%>%
  mutate(.,
         new_word = factor(word, levels=word[order(n,decreasing=F)]))


myplotsIL <- lapply(split(clean_text_IL,list(clean_text_IL$country,clean_text_IL$il_condition)), function(x){
  x$word <- factor(x$word, levels=x$word[order(x$n,decreasing=F)])
  x<-x%>%
    mutate(.,
           il_condition=case_when(
             il_condition=="no_law"~"No Law",
             il_condition=="law"~"Law"))
  x$country_lab<-ifelse(x$il_condition=="Law",x$country,"")
  
  #make the plot
  p <- ggplot(x, aes(x = word, y = n, width=0.75)) +
    geom_bar(stat = "identity") + #to force all levels to be considered
    ylim(0,400)+
    theme_bw()+
    theme(legend.position="none")+
      labs(y="count",x=unique(x$country_lab),title = unique(x$il_condition))+
    coord_flip()
})

myplotsIL <- myplotsIL[order(names(myplotsIL))]
a<-do.call(arrangeGrob,(c(myplotsIL, ncol=2)))
ggsave("figures/appendix/words_IL.pdf",a, width = 8, height = 14)

#IPE
clean_text_IPE <- all_countries_clean %>%
  filter(info_leak_ipe_1_TEXT!="") %>%
  dplyr::select(info_leak_ipe_1_TEXT,country,ipe_condition)%>%
  unnest_tokens(word, info_leak_ipe_1_TEXT)%>% 
  anti_join(stop_words, by = "word") %>%
  group_by(country,ipe_condition) %>%
  count(word) %>% 
  top_n(5, n) %>%
  arrange(desc(n), .by_group = TRUE)%>%
  mutate(.,
         new_word = factor(word, levels=word[order(n,decreasing=F)]))


myplotsIPE <- lapply(split(clean_text_IPE,list(clean_text_IPE$country,clean_text_IPE$ipe_condition)), function(x){
  x$word <- factor(x$word, levels=x$word[order(x$n,decreasing=F)])
  x<-x%>%
    mutate(.,
           ipe_condition=case_when(
             ipe_condition=="easier"~"Easier",
             ipe_condition=="harder"~"Harder"))
  x$country_lab<-ifelse(x$ipe_condition=="Easier",x$country,"")
  
  #make the plot
  p <- ggplot(x, aes(x = word, y = n, width=0.75)) +
    geom_bar(stat = "identity") + #to force all levels to be considered
    ylim(0,400)+
    theme_bw()+
    theme(legend.position="none")+
    labs(y="count", x=unique(x$country_lab),
         title = unique(x$ipe_condition))+
    coord_flip()
})

myplotsIPE <- myplotsIPE[order(names(myplotsIPE))]

a<-do.call(arrangeGrob,(c(myplotsIPE, ncol=2)))
ggsave("figures/appendix/words_IPE.pdf",a, width = 8, height = 14)


#Plot means of outcomes by conditions by country
##Estimate Country Specific means by condition:
conditions_DP_bycountry <-
  all_countries_clean %>%
  mutate(.,Conditions=case_when(dem_condition=="democracy"~"Treatment",
                                   dem_condition=="non_democracy"~"Control")) %>% 
  filter(Conditions!="")%>%
  group_by(country,Conditions) %>%
  do(tidy(lm_robust(outcome1_dem_new ~ 1, data = .)))%>%
  mutate(.,
         meta = "Country Specific Replication",
         study = "Democratic Peace \n (Support Attack)")

##Estimate for all countries
conditions_DP_all <- 
  all_countries_clean %>% 
  mutate(.,Conditions=case_when(dem_condition=="democracy"~"Treatment",
                                dem_condition=="non_democracy"~"Control")) %>% 
  filter(Conditions!="")%>%
  group_by(Conditions) %>%
  do(tidy(lm_robust(outcome1_dem_new ~ 1, data = .)))%>%
  mutate(.,
         country = "All countries",
         meta = "All",
         study = "Democratic Peace \n (Support Attack)")

##Estimate original means by condition:
conditions_DP_orig <- 
  dem_peace_orig %>%
  mutate(.,Conditions=case_when(democ==1~"Treatment",
                           democ==0~"Control")) %>%
  group_by(Conditions) %>%
  do(tidy(lm_robust(strike ~ 1, data = .)))%>%
  mutate(.,
         meta = "Original",
         country = "USA",
         study = "Democratic Peace \n (Support Attack)")

##Audience costs
##Country specific means

##Estimate Country Specific means by condition:

conditions_AC_bycountry <-
  all_countries_clean %>%
  mutate(.,Conditions=case_when(aud_condition=="not_engage"~"Treatment",
                                   aud_condition=="stay_out"~"Control")) %>% 
  filter(Conditions!="")%>%
  group_by(country,Conditions) %>%
  do(tidy(lm_robust(outcome_aud_new ~ 1, data = .)))%>%
  mutate(.,
         meta = "Country Specific Replication",
         study = "Audience Costs \n (Leader Approval)")

##Estimate for all countries
conditions_AC_all <- 
  all_countries_clean %>% 
  mutate(.,Conditions=case_when(aud_condition=="not_engage"~"Treatment",
                                   aud_condition=="stay_out"~"Control")) %>% 
  filter(Conditions!="")%>%
  group_by(Conditions) %>%
  do(tidy(lm_robust(outcome_aud_new ~ 1, data = .)))%>%
  mutate(.,
         country = "All countries",
         meta = "All",
         study = "Audience Costs \n (Leader Approval)")

##Estimate original means by condition:
conditions_AC_orig <- 
  audience_orig %>%
  mutate(.,
         Conditions = case_when(
           Strategy == 3 ~ "Treatment",
           Strategy == 1 ~ "Control"),
         outcome = case_when(
           totApprov1== -3 ~1,
           totApprov1== -2 ~2,
           totApprov1== -1 ~3,
           totApprov1== 0~4,
           totApprov1== 1~5,
           totApprov1== 2~6,
           totApprov1== 3~7)) %>%
  group_by(Conditions)%>%
  do(tidy(lm_robust(outcome~1, data = .))) %>% 
  mutate(.,
         meta = "Original",
         country = "USA",
         study = "Audience Costs \n (Leader Approval)")%>%
  drop_na()


##International Law
##Country specific means

##Estimate Country Specific means by condition:
conditions_IL_bycountry <-
  all_countries_clean %>%
  mutate(.,Conditions=case_when(il_condition=="law"~"Treatment",
                                   il_condition=="no_law"~"Control")) %>% 
  filter(Conditions!="")%>%
  group_by(country,Conditions) %>%
  do(tidy(lm_robust(outcome_il_new ~ 1, data = .)))%>%
  mutate(.,
         meta = "Country Specific Replication",
         study = "International Law \n (Support Torture)")

##Estimate for all
conditions_IL_all <-
  all_countries_clean %>%
  mutate(.,Conditions=case_when(il_condition=="law"~"Treatment",
                                  il_condition=="no_law"~"Control")) %>% 
  filter(Conditions!="")%>%
  group_by(Conditions) %>%
  do(tidy(lm_robust(outcome_il_new ~ 1, data = .)))%>%
  mutate(.,
         country = "All countries",
         meta = "All",
         study = "International Law \n (Support Torture)")

##Estimate original:

conditions_IL_orig <-
  il_orig %>% 
  mutate(.,Conditions=case_when(treaty==1~"Treatment",
                                  treaty==0~"Control")) %>%
  group_by(Conditions)%>%
  do(tidy(lm_robust(torture7~1,data=.))) %>% 
  mutate(.,
         meta = "Original",
         country = "USA",
         study = "International Law \n (Support Torture)")


##IPE
##Country specific means
##Estimate Country Specific means by condition:
conditions_IPE_bycountry <-
  all_countries_clean %>%
  mutate(.,Conditions=case_when(ipe_condition=="harder"~"Treatment",
                                  ipe_condition=="easier"~"Control")) %>% 
  filter(Conditions!="")%>%
  group_by(country,Conditions) %>%
  do(tidy(lm_robust(outcome_ipe_new ~ 1, data = .)))%>%
  mutate(.,
         meta = "Country Specific Replication",
         study = "Reciprocity FDI \n (Support harder barriers)")

##Estimate for all
conditions_IPE_all<-
  all_countries_clean %>%
  mutate(.,Conditions=case_when(ipe_condition=="harder"~"Treatment",
                                   ipe_condition=="easier"~"Control")) %>% 
  filter(Conditions!="")%>%
  group_by(Conditions) %>%
  do(tidy(lm_robust(outcome_ipe_new ~ 1, data = .)))%>%
  mutate(.,
         country = "All countries",
         meta = "All",
         study = "Reciprocity FDI \n (Support harder barriers)")

##Estimate Original
conditions_IPE_orig<-
  reciprocity_orig%>%
  mutate(.,Conditions= case_when(current==0~"Control",
                                    current==6~"Treatment"),
         outcome = case_when(pref==0.00~1,
                             pref==0.25~2,
                             pref==0.50~3,
                             pref==0.75~4,
                             pref==1.00~5))%>%
  group_by(Conditions) %>%
  drop_na(outcome)%>%
  do(tidy(lm_robust(outcome ~ 1, data = .)))%>%
  mutate(.,
         meta = "Original",
         country = "USA",
         study = "Reciprocity FDI \n (Support harder barriers)")


##combine all means for all experiments

conditions_df <-
  bind_rows(conditions_DP_bycountry,conditions_DP_all,conditions_DP_orig,
            conditions_AC_bycountry,conditions_AC_all,conditions_AC_orig,
            conditions_IL_bycountry,conditions_IL_all,conditions_IL_orig,
            conditions_IPE_bycountry,conditions_IPE_all,conditions_IPE_orig)

conditions_df$country <- 
  factor(conditions_df$country,
         levels = c("USA","Nigeria","Japan", "India",
                    "Israel","Germany", "Brazil", "All countries"))

conditions_df$Conditions <- 
  factor(conditions_df$Conditions,
         levels = c("Treatment","Control"))

conditions_df$meta <- 
  factor(conditions_df$meta,
         levels = c("All","Country Specific Replication","Original"))

ggplot(conditions_df, aes(x = as.numeric(estimate), y =as.factor(country), 
                               color =Conditions, shape = Conditions)) +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high), size = 0.3) +
  labs(x = "Mean",
       y = "") +
  facet_grid(cols = vars(study), rows = vars(meta),
             scales="free_y", space = "free")+
  scale_color_manual(values = c("black","grey60"))+
  theme_bw()+
  theme(strip.background =element_rect(fill="dodgerblue4"),
        strip.text = element_text(colour = 'white'),
        panel.grid.major = element_blank(),
        #legend.position = "none",
        strip.text.x = element_text(size = 12),
        strip.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 12),
        axis.text.x = element_text(size = 11),
        axis.text.y = element_text(size = 11),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
       strip.placement='outside')
ggsave("figures/appendix/analysis_byconditions.pdf", width = 11, height = 6)


#External validity bias----
#change covariates to numeric

covariates <- c("gender_fe", "age", "dem_norm", "hawk", "legal", 
            "fear_1", "fear_2","voting_fe", "ideology", "education_fe") ##set covariates

#create plot per country per experiment
#Brazil
bra_DP<-all_countries_clean%>%
  filter(country=="Brazil")%>%
  drop_na(outcome1_dem_scaled)%>%
  exr(outcome = "outcome1_dem_scaled", 
      treatment = "DP_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

summary(bra_DP)

bra_AC<-all_countries_clean%>%
  filter(country=="Brazil")%>%
  drop_na(outcome_aud_scaled)%>%
  exr(outcome = "outcome_aud_scaled", 
      treatment = "AC_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)
summary(bra_AC)

bra_IL<-all_countries_clean%>%
  filter(country=="Brazil")%>%
  drop_na(outcome_il_scaled)%>%
  exr(outcome = "outcome_il_scaled", 
      treatment = "IL_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)
summary(bra_IL)

bra_IPE<-all_countries_clean%>%
  filter(country=="Brazil")%>%
  drop_na(outcome_ipe_scaled)%>%
  exr(outcome = "outcome_ipe_scaled", 
      treatment = "IPE_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)
summary(bra_IPE)

#Germany
ger_DP<-all_countries_clean%>%
  filter(country=="Germany")%>%
  drop_na(outcome1_dem_scaled)%>%
  exr(outcome = "outcome1_dem_scaled", 
      treatment = "DP_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)
summary(ger_DP)


ger_AC<-all_countries_clean%>%
  filter(country=="Germany")%>%
  drop_na(outcome_aud_scaled)%>%
  exr(outcome = "outcome_aud_scaled", 
      treatment = "AC_treat", 
      covariates = covariates,
      data = .,
      sate_estimate = NULL)


ger_IL<-all_countries_clean%>%
  filter(country=="Germany")%>%
  drop_na(outcome_il_scaled)%>%
  exr(outcome = "outcome_il_scaled", 
      treatment = "IL_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

ger_IPE<-all_countries_clean%>%
  filter(country=="Germany")%>%
  drop_na(outcome_ipe_scaled)%>%
  exr(outcome = "outcome_ipe_scaled", 
      treatment = "IPE_treat", 
      covariates =  covariates,
      data = .,
      sate_estimate = NULL)

#India
ind_DP<-all_countries_clean%>%
  filter(country=="India")%>%
  drop_na(outcome1_dem_scaled)%>%
  exr(outcome = "outcome1_dem_scaled", 
      treatment = "DP_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

ind_AC<-all_countries_clean%>%
  filter(country=="India")%>%
  drop_na(outcome_aud_scaled)%>%
  exr(outcome = "outcome_aud_scaled", 
      treatment = "AC_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

ind_IL<-all_countries_clean%>%
  filter(country=="India")%>%
  drop_na(outcome_il_scaled)%>%
  exr(outcome = "outcome_il_scaled", 
      treatment = "IL_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

ind_IPE<-all_countries_clean%>%
  filter(country=="India")%>%
  drop_na(outcome_ipe_scaled)%>%
  exr(outcome = "outcome_ipe_scaled", 
      treatment = "IPE_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

#Israel
isl_DP<-all_countries_clean%>%
  filter(country=="Israel")%>%
  drop_na(outcome1_dem_scaled)%>%
  exr(outcome = "outcome1_dem_scaled", 
      treatment = "DP_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

isl_AC<-all_countries_clean%>%
  filter(country=="Israel")%>%
  drop_na(outcome_aud_scaled)%>%
  exr(outcome = "outcome_aud_scaled", 
      treatment = "AC_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

isl_IL<-all_countries_clean%>%
  filter(country=="Israel")%>%
  drop_na(outcome_il_scaled)%>%
  exr(outcome = "outcome_il_scaled", 
      treatment = "IL_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)


#remove fear_2 to increase optimization
isl_IPE<-all_countries_clean%>%
  filter(country=="Israel")%>%
  drop_na(outcome_ipe_scaled)%>%
  exr(outcome = "outcome_ipe_scaled", 
      treatment = "IPE_treat", 
      covariates = c("gender_fe", "age", "dem_norm", "hawk", "legal", 
                     "voting_fe", "fear_1", "ideology", "education_fe"),
      data = .,
      sate_estimate = NULL)


#Japan
jap_DP<-all_countries_clean%>%
  filter(country=="Japan")%>%
  drop_na(outcome1_dem_scaled)%>%
  exr(outcome = "outcome1_dem_scaled", 
      treatment = "DP_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

jap_AC<-all_countries_clean%>%
  filter(country=="Japan")%>%
  drop_na(outcome_aud_scaled)%>%
  exr(outcome = "outcome_aud_scaled", 
      treatment = "AC_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

jap_IL<-all_countries_clean%>%
  filter(country=="Japan")%>%
  drop_na(outcome_il_scaled)%>%
  exr(outcome = "outcome_il_scaled", 
      treatment = "IL_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

jap_IPE<-all_countries_clean%>%
  filter(country=="Japan")%>%
  drop_na(outcome_ipe_scaled)%>%
  exr(outcome = "outcome_ipe_scaled", 
      treatment = "IPE_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

#Nigeria
nig_DP<-all_countries_clean%>%
  filter(country=="Nigeria")%>%
  drop_na(outcome1_dem_scaled)%>%
  exr(outcome = "outcome1_dem_scaled", 
      treatment = "DP_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

nig_AC<-all_countries_clean%>%
  filter(country=="Nigeria")%>%
  drop_na(outcome_aud_scaled)%>%
  exr(outcome = "outcome_aud_scaled", 
      treatment = "AC_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

nig_IL<-all_countries_clean%>%
  filter(country=="Nigeria")%>%
  drop_na(outcome_il_scaled)%>%
  exr(outcome = "outcome_il_scaled", 
      treatment = "IL_treat", 
      covariates =  covariates, 
      data = .,
      sate_estimate = NULL)

nig_IPE<-all_countries_clean%>%
  filter(country=="Nigeria")%>%
  drop_na(outcome_ipe_scaled)%>%
  exr(outcome = "outcome_ipe_scaled", 
      treatment = "IPE_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

#USA
us_DP<-all_countries_clean%>%
  filter(country=="USA")%>%
  drop_na(outcome1_dem_scaled)%>%
  exr(outcome = "outcome1_dem_scaled", 
      treatment = "DP_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

us_AC<-all_countries_clean%>%
  filter(country=="USA")%>%
  drop_na(outcome_aud_scaled)%>%
  exr(outcome = "outcome_aud_scaled", 
      treatment = "AC_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

us_IL<-all_countries_clean%>%
  filter(country=="USA")%>%
  drop_na(outcome_il_scaled)%>%
  exr(outcome = "outcome_il_scaled", 
      treatment = "IL_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)

us_IPE<-all_countries_clean%>%
  filter(country=="USA")%>%
  drop_na(outcome_ipe_scaled)%>%
  exr(outcome = "outcome_ipe_scaled", 
      treatment = "IPE_treat", 
      covariates = covariates, 
      data = .,
      sate_estimate = NULL)


pdf(file="figures/appendix/external_val.pdf", width = 15, height = 17)
par(mfrow=c(7,4))

plot(bra_AC)
mtext("Audience Costs (Brazil)",side=3)
plot(bra_DP)
mtext("Democratic Peace (Brazil)",side=3)
plot(bra_IL)
mtext("International Law (Brazil)",side=3)
plot(bra_IPE)
mtext("FDI Reciprovity (Brazil)",side=3)

plot(ger_AC)
mtext("Audience Costs (Germany)",side=3)
plot(ger_DP)
mtext("Democratic Peace (Germany)",side=3)
plot(ger_IL)
box(col = "red")  
mtext("International Law (Germany)",side=3)
plot(ger_IPE)
mtext("FDI Reciprovity (Germany)",side=3)

plot(ind_AC)
mtext("Audience Costs (India)",side=3)
plot(ind_DP)
box(col = "red")  
mtext("Democratic Peace (India)",side=3)
plot(ind_IL)
box(col = "red")  
mtext("International Law (India)",side=3)
plot(ind_IPE)
mtext("FDI Reciprovity (India)",side=3)

plot(isl_AC)
mtext("Audience Costs (Israel)",side=3)
plot(isl_DP)
mtext("Democratic Peace (Israel)",side=3)
plot(isl_IL)
box(col = "red")  
mtext("International Law (Israel)",side=3)
plot(isl_IPE)
mtext("FDI Reciprovity (Israel)",side=3)

plot(jap_AC)
mtext("Audience Costs (Japan)",side=3)
plot(jap_DP)
box(col = "red")  
mtext("Democratic Peace (Japan)",side=3)
plot(jap_IL)
mtext("International Law (Japan)",side=3)
plot(jap_IPE)
box(col = "red")  
mtext("FDI Reciprovity (Japan)",side=3)

plot(nig_AC)
mtext("Audience Costs (Nigeria)",side=3)
plot(nig_DP)
box(col = "red")  
mtext("Democratic Peace (Nigeria)",side=3)
plot(nig_IL)
mtext("International Law (Nigeria)",side=3)
plot(nig_IPE)
mtext("FDI Reciprovity (Nigeria)",side=3)

plot(us_AC)
mtext("Audience Costs (USA)",side=3)
plot(us_DP)
mtext("Democratic Peace (USA)",side=3)
plot(us_IL)
mtext("International Law (USA)",side=3)
plot(us_IPE)
mtext("FDI Reciprovity (USA)",side=3)
dev.off()


#Probing the null in India (DP) ----

#remove subjects who failed manipulations of other studies, then compare to India DP
#AC manipulation
all_countries_clean$AC_man<- ifelse(all_countries_clean$manipulation_aud=="1"&
                                      all_countries_clean$aud_condition == "stay_out",1,NA)
all_countries_clean$AC_man<- ifelse(all_countries_clean$manipulation_aud=="2"&
                                      all_countries_clean$aud_condition == "not_engage",1,
                                    all_countries_clean$AC_man)
all_countries_clean$AC_man<- ifelse(all_countries_clean$manipulation_aud=="3"&
                                      all_countries_clean$aud_condition == "engage",1,
                                    all_countries_clean$AC_man)

#IL manipulation
all_countries_clean$IL_man<- ifelse(all_countries_clean$manipulation_il=="1"&
                                      all_countries_clean$il_condition=="law",1,NA)
all_countries_clean$IL_man<- ifelse(all_countries_clean$manipulation_il=="2"&
                                      all_countries_clean$il_condition=="no_law",1,
                                    all_countries_clean$IL_man)

#IPE manipulation
all_countries_clean$IPE_man<- ifelse(all_countries_clean$manipulation_ipe=="1"&
                                       all_countries_clean$ipe_condition=="easier",1,NA)
all_countries_clean$IPE_man<- ifelse(all_countries_clean$manipulation_ipe=="2"&
                                       all_countries_clean$ipe_condition=="harder",1,
                                     all_countries_clean$IPE_man)

india_attentive<-all_countries_clean%>%
  filter(country=="India"&
           AC_man==1&
           IL_man==1&
           IPE_man==1)

india_attentive<-lm(outcome1_dem_scaled~DP_treat,india_attentive)

texreg(india_attentive,
       label = "tab:india_attentive",
       caption.above = T,
       no.margin = T,
       include.ci = FALSE,
       custom.header = list("Support attack" = 1),
       custom.coef.map = list("DP_treat" = "Democracy"),
       stars = 0.05,
       digits = 3,
       include.rsquared = F, include.rmse = F,
       caption = "Screening out failed manipulation from other studies (India DP)",
       font.size = "tiny",
       table = T,
       custom.note = "%stars.",
       file = "tables/india_attentive.tex")

